UNCLASSIFIED 


SECURITY  CLASSIRCATION  OF  THIS  PAGE 


REPORT  DOCUMENTATION  PAGE 


Form  Appf0¥0d 
0MB  Mo.  0704^13$ 


la.  REPORT  SECURITY  CLASSIRCATON 

UnclassiOed 


lb.  RESTRICTIVE  MARKINGS 


2a.  SECURITY  CLASSIFICATION  AUTHORITY 


3.  DISTRIBUTION/AVAIIABILITY  OF  REPORT 


2b.  OECLASSIFICATIOf^CXWVNGRAOING  SCHEDULE 


Approved  for  public  release;  distribution  unlimited. 


4.  PERFORMING  ORGANIZATION  REPORT  NUMBER(S) 

HDL-TR-2184 


S.  MONITORING  ORGANIZATION  REPORT  NUMeER(S) 


6a.  NAME  OF  PERFORMING  ORGANIZATON 

Harry  Diamond  Laboratories 


6b.  OFFICE  SYluBOL 
(H  appttcablo) 

SLCHD-NW-EP 


7a.  NAME  OF  MONITORING  ORGANIZATION 


6c  ADDRESS  (C»/,  Stale,  and  ZIP  Code) 

2800  Powder  Mill  Road 
Adelphi.MD  20783-1 197 


7b.  ADDRESS  (City.  State,  and  ZIP  Code) 


Sa.  NAME  OF  FUNOING/SPONSORING 
ORGANIZATION 

U.S.  Army  Laboratory  Command 


7b.  OFFICE  SYMBOL 
^applicable) 

AMSLC 


9.  PROCUREMENT  INSTRUMENT  IDENTIFICATION  NUMBER 


6c.  ADDRESS  (Cly.  State,  and  ZIP  Code) 

2800  Powder  Mill  Road 
Adelphi.MD  20783-1 145 


10.  SOURCE  OF  FUNI 


PROGRAM 
ELEMENT  NO. 

6.11.02 


DING  NUMBERS _ I 

PROJECT 

NO. 

TASK 

NO. 

WORK  UNIT 
ACCESSION  NO. 

11.  TITLE  (Include  Security  Classiication) 

Mobility  of  Nonequilibrium  Conduction  Electrons  in  Air 


12.  PERSONAL  AUTHOR(S) 

William  T.  Wyatt,  Jr.,  and  Christopher  S.  Kenyon 


13a.  TYPE  OF  REPORT 

Final 


13b.  TIME  COVERED 

FROM  Oct  88  TO  Sept  89 


14.  DATE  OF  REPORT  (Year,  Month,  Day) 

June  1990 


15.  PAGE  COUNT 

46 


16.  SUPPLEMENTARY  NOTATION 

AMS  code:  61 1 102.H44001 1 


7 


1  17.  COSATI  COOES  (  /  1 

RELD 

GROUP 

SUB-GRObP — ^7 

20 

9 

18.  SUBJECT  TEFIMS  (Continue  on  m/eree  H  neceaaaty  and  tdentty  by  btock  number) 

Plasma,  conductivity,  transients,  variational  principle,  electron , 


\l 


19.  ABSTRACT  (CantinuB  on  mvofoe  0  nocessaiy  and  kMntiy  by  block  numbor) 


An  integro-differential  equation  is  derived  describing  the  time  evolution  of  the  electron  energy  distribution  function  in  a 
gas  in  a  transient  electric  field.  A  finite-difference  approximation  and  a  time-stepping  algorithm  using  matrix-vector 
techniques  are  devised  to  solve  the  integro-differential  equation.  An  extensive  published  cross  section  set  is  added  to  the 
time-stepping  algorithm,  enabling  one  to  calculate  realistic  energy  distribution  functions  in  air.  Also  described  is  a 
projection  method  that  allows  computation  of  electron  swarm  properties  such  as'*'collision  volume^nversely  related  to 
mobility)  using  measured  data  for  equilibrium  distributions.  This  projection  method  is  shown  to  have  desirable  variational 
properties  in  calculating  the  swarm  properties.  W^t^ompare  orthogonal  and  nonnegative  projection  methods  with 
conventional  techniques  for  their  ability  to  reduce  error  iif  estimates  of  collision  volume  arisipg  from  cross-section-induced 
error  in  the  energy  distribution  function.  We  found  that\single-vector  projection  with  ^’d^up-'^bf  approximation  error 
performs  best,  leading  to  significant  reduction  of  collision' volume  error  when  the  distribution  function  maximum  occurs 


20.  DISTRiaUTIONIAVAILABILnYOF  ABSTRACT 
n  UNCLASSIFIEDAJNLIMRED  PsAMEASRPT.  □  DTIC  USERS 


21.  ABSTRACT  SECURITY  CLASSIFICATION 

Unclassified 


22a.  NAME  OF  RESPONSIBLE  INDIVIDUAL 

William  T.  Wyatt,  Jr. 


22b.  TELEPHONE  (Include  Area  Code) 

(703)  490-2303 


22e.  OFFICE  SYMBOL 

SLCHD-NW-EP 


00  Form  1473,  JUN  86 


Prerioue  edUom  ate  obaoleia. 


SECURITY  CLASSIFICATION  OF  THIS  PAGE 

UNCLASSIFIED 


UNCLASSIFIED 


SECURITY  CU3SIFCATION  OF  THIS  PAGE 


19.  Abstract  cont’d 

at  not  more  than  about  0 . 1  e  V .  Other  projections  perform  poorly  compared  to  conventional  methods  of  calculating  swarm 
collision  volume,  because  of  inherent  dependence  of  the  projection  basis  of  equilibrium  distributions  and  resulting  ill- 
conditioning  of  the  projection  matrix.  /  /  ‘  // 


AcCcSK’'  fO’ 

NflS  CRAAI 
DUO  lAb 

UnannOit'Ccd 
Justif.cat-O'' _ 

By _ 

Distribution/ 


/.Tilabiiity  Codes 

Ava'i  and /or 

Dist 

Special 

lI _ ! 

UNCLASSIFIED 


Table  of  Contents 


1  Introduction  and  Motivation  . 1 

1.1  Nonequilibrium  .Mobility  by  Projection  onto  an  Pquilibrium  Basis  . I 

1.2  Relationship  to  a  Variational  Principle  . 3 

2  Solution  Method  for  F.ncrgy  Distribution  Function  . 5 

2.1  Fnergy  Diffusion  1-quation  . 5 

2.2  Finite  Mesh  Model  of  the  Fnergy  Diffusion  liquation  . 12 

2.3  Flastic  Scattering  . 14 

2.4  Unequally  Spaced  Fnergy  Mesh  . 17 

2.5  Fxact  Solutions  and  Equilibrium  Finite  Mesh  Solutions  Compared  ....  18 

2.6  Inelastic  Scattering  . 23 

2.7  Equilibrium  Energy  Distribution  Compared  to  Previous  Results  . 25 

2.8  Equilibrium  Transport  Properties  Compared  to  Vleasured  Data  . 26 

3  Analysis  of  l-quilibrium  Projeetion  Method  . 27 

3.1  Developing  the  Basis  of  Equilibrium  Energy  Distributions  . 27 

3.2  Obtaining  Time-Dependent  Energy  Distributions  . 33 

3.2  Projecting  Fi me- Dependent  Distributions  onto  an  Equilibrium  Basis  ...  34 

3.3  Defining  a  Test  of  the  Equilibrium  Projection  .Vlethod  . 35 

3.4  Orthogonal  Projection  Results  for  a  Time- Dependent  Distribution  ....  37 

3.5  Nonnegative  Projection  Method  . 37 

4  Discussion  . 40 

5  Conclusions  . 42 

6  Recommendations  . 43 

References  . 44 

Distribution  . 45 

List  of  Illustrations 

Figure  1.  Finite  mesh  equilibrium  results  for  Maxwellian  distribution  . 21 

I'igurc  2.  Finite  mesh  equilibrium  results  for  Druyvesteyn  distribution  . 22 

I'igurc  3.  Finite  mesh  results  and  first  two  I^gcndre  coefficients  compared  ...  25 

Figure  4.  Equilibrium  energy  distributions  as  trial  basis  vectors  . 30 

Figure  5.  Orthonormal  basis  vectors  derived  from  equilibrium  energy  distrib¬ 
utions  . 31 

Figure  6.  Nonequilibrium  energy  distributions  for  step-function  electric  field  .  .  32 

Figure  7.  Nonequilibrium  energy’  distributions  for  square-wave  electric  field  ...  33 

Figure  8.  Spectral  projections  for  step-function  electric  field  . 35 

Figure  9.  Spectral  projections  for  square-wave  electric  field  . 36 

Figure  10.  f  'rror  in  calculated  ensemble  collision  volume  for  step-function  E-ficld  39 

Figure  1 1.  Error  in  calculated  ensemble  collision  volume  for  square-wave  F  field  40 

List  of  Tables 

Table  1.  Calculated  and  Measured  Swarm  Parameters  in  Nitrogen  Compared  .  26 


1  Introduction  and  Motivation 


1.1  Nonequilibrium  Mobility  by  Projection  onto  an  Equilibrium  Basis 


The  mobility  of  conduction  electrons  in  air  subjected  to  an  electric  field 
changing  rapidly  in  a  nanosecond  or  less  is  difficult  to  measure  or  calculate. 
The  conduction  electrons  do  not  attain  an  equilibrium  state  in  such  a  short 
time.  Steady-state  “swarm”  measurements  are  readily  available,  but  only 
appro.ximate  nonequilibrium  values.  Approaches  using  energy-dependent 
momentum  exchange  cross-section  data  with  a  conduction  electron  velocity 
distribution  function  arc  hampered  by  incomplete  cross-section  data.  In  this 
work,  we  project  the  time-dependent  energy  distribution  function  onto  a 
basis  composed  of  equilibrium  energy  distribution  functions  and  use  equi¬ 
librium  mobilities  to  construct  a  nonequilibrium  mobility.  The  method  is 
evaluated  as  a  possible  means  of  avoiding  error  in  the  energy  distribution 
function  arising  from  incomplete  or  erroneous  cross-section  data. 

Conduction  electron  mobility  can  be  expressed  in  terms  of  an  effective  mo¬ 
mentum  exchange  cross-section  as'  [pp  72,  115,  206,  527] 


qW  10'" 

(\a) 

9 

muNQ^  ’ 

(Ih) 

where 

fx  =  electron  mobility  (in  m^fV  s), 

Q„  —  effective  momentum  exchange  cross-section  (in  m^molcculc), 
q  =  magnitude  of  the  electron  charge  (1.60  X  10  ”  C), 
iV  =  gas  molecular  weight  (28.8  kg/Molc,  1  Mole  =  1000  mole), 
m  =  electron  mass  (9.1 1  x  10  ^' kg), 

u  =  electron  speed  (m/s)  for  electron  kinetic  energy  E  (in  J), 

A„  =  Avogadro's  number  (6.03  x  lO^^f^olc), 

p  =  gas  density  (in  kg/m^),  and 

iV  =  gas  molecular  density  (in  molcculcs/m’). 

The  macroscopic  conduction  electron  mobility,  averaged  over  the  ensemble 
electron  energies  (which  we  denote  by  a  caret),  is  written  as 


_ 9 _ 

mNju{E:)Q^{E)f[E)dE 


(2) 


'  Huxley,  L.  G.  M..  and  R.  VV.  Crompton,  The  Diffusion  and  Drift  of  Electrons  in  Gases,  John 
Wiley  and  Sons,  New  York  (1974). 
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where  J{E)  is  the  time-independent  (equilibrium)  energy  distribution  of 
electron  kinetic  energies.  We  may,  for  convenience,  parameterize  as  a 

function  of  the  ensemble  average  electron  energy  /s,  which  is  a  monotonic 
function  of  the  ambient  (time-independent)  electric  field  strength.  Let  us 
define 

K{E)  =  u{i:)Q„{E)  (3) 

A  A 

and  the  ensemble  average  as 

k{iC)  =  ju{i-:)Q^{iW-:)di-:.  (4) 

The  variables  K  and  K,  which  we  shall  call  “collision  volume,”  arc  in  units 
of  volume  per  unit  time  and,  when  multiplied  by  the  density  of  target  mole¬ 
cules,  yield  the  momentum  exchange  collision  frequency  v„. 

The  ensemble  average  collision  volume  K  can  be  written  as  a  function  of 
time 


m = J  dE  (5) 

where  the  integral  is  carried  over  all  electron  energies,  and  _/(/'./)  is  the  time- 
dependent  energy  distribution  function  (in  electrons  per  unit  energy).  We 
assert  without  proof  that  an  infinite  basis  exists  complete  for  any  J{E,t), 
consisting  of  the  family  of  equilibrium  energy  distribution  functions  (/>(£,  A) 
for  all  nonnegative  values  of  the  parameter  A.  The  parameter  A  may  be 
chosen  to  be  some  monotonic  characteristic  of  the  distribution,  such  as 
magnitude  of  the  electric  field  strength  or  average  energy.  We  discretize  the 
basis  for  parameter  values  A  =  A,  and  approximate  by 

y(/:,O  =  Xa,(/)0(£.A,),  (6) 

i 

which  defines  time-dependent  cocfilcicnts  a,{t).  (  This  definition  docs  not 
uniquely  define  the  a,(/),  however,  because  the  basis  is  not  orthogonal. 
Orthogonality  will  be  considered  in  section  3.1.)  Substituting  this  basis  into 
equation  (5)  we  obtain 

ki)  =  K{I-:)a,{i)4>{E,  A,)  dE, 

*  (7) 

=  ^a,(/)j  K{imE,\;)dE. 

i 

The  ensemble  average  collision  volume  for  any  equilibrium  distribution 
having  parameter  value  A,  is 

AW  =  j  A(A)0(£,A,)rfA.  (8) 

The  ensemble  average  collision  volume  K  can  be  easily  obtained  from 
“swarm”  measurements  of  electron  mobility,  which  are  aimed  at  taking  data 
for  equilibrium  energy  distributions.  By  substitution  of  equation  (8)  into 
equation  (7)  we  obtain 


2 


(9) 


i 

Thus  the  time-dependent  ensemble  average  collision  volume  can  be  obtained 
from  equilibrium  ensemble  average  collision  volumes  and  a  knowledge  of  the 
a,(r)  describing  the  time-dependent  energy  distribution.  Alternatively,  the 
same  data  can  be  derived  from  the  (energy-dependent)  momentum  exchange 
collision  cros.s-scction  and  the  time-dependent  energy  distribution  itself 
Unfortunately  these  quantities  arc  not  known  accurately  for  all  energies. 

The  motivation  for  the  current  work  is  the  possibility  that  the  coefficients 
«,(/)  can  be  determined  more  accurately  than  either  ]{E,i)  or  A,),  and 
also  more  accurately  than  the  collision  volume  data  K{1?\  upon  which  they 
arc  based.  Under  certain  conditions,  the  cocfTicicnts  a,{i)  arc  stationary  with 
respect  to  variation  of  K{1^  (due  to  a  variational  principle),  for  example, 
this  occurs  when  the  time-dependent  electron  energy  distribution  function  is 
close  to  an  equilibrium  energy  distribution  function.  More  generally,  sta¬ 
tionary  behavior  can  be  expected  when  the  time-dependent  function  can  be 
approximated  by  a  set  of  linearly  independent  equilibrium  functions.  If  in¬ 
deed  the  a{t)  can  be  determined  more  accurately,  then  the  ensemble  average 

collision  volume  A'(A,)  and  the  ensemble  average  mobility  A  can  also  be  de¬ 
termined  more  accurately  than  by  direct  integration  of  equation  (5). 

,2  Relationship  to  a  Variational  Principle 

Define  a  linear  operator  L  that  advances  the  energy  distribution  function 
J[E,t)  one  time-step: 

Ly(£,/)  =y(/:,/ +  A/).  (lo) 

The  equilibrium  energy  distribution  function  is  defined  by 

Lg(/r)  =  g(E)  (11) 

for  a  fixed  electric  field  strength.  Clearly  giE)  is  an  eigenvector  of  L  associ¬ 
ated  with  the  eigenvalue  unity.  A  well-known  result^  shows  that  eigenvalues 
may  be  derived  from  a  variational  principle  as 

f  [Ly(x)]/(x)^x 

^=^—r -  (>2) 

J  f\x)  dx 

and  that  this  expression  for  A  is  stationary  (second-order  accuracy)  with  re¬ 
spect  to  variations  inj[x). 

Identify  the  g(f'.')  that  has  the  largest  projection  onto  i.  e.,  with  the 

largest  «(/)  defined  by 


AE,i)g{E)dE 

=  ^ - .  (13) 

g\l^dE 


^ Morse,  P.  M.,  and  H.  Feshbach,  Methods  of  Theoretical  Physics,  McGraw-Hill  Book 
Company,  Inc.,  New  York  (1953).  Volume  2,  pp  1 108f. 
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Repeated  applications  of  L  lo  J{E,t)  produce  ncwJ{E,t  +  At)  that  grow  closer 
to  the  equilibrium  solution  g{E) 

L"j{E,t)  =  git’)  +  Sg{E,l),  (14) 

where  5g{E,t)  is  a  differential  variation  of  function  g,  such  that 

lim  =  0.  (15) 

n-*oo 

Now  define  an  operator  L  '  that  is  the  inverse  of  L.  We  observe  that  re¬ 
peated  applications  of  L  '  to  g(I:)  -I-  Sg{E,t)  produce 

(L-')''[g(£)  -f  SgiE,t)^  =y(£,/).  (16) 

If  we  now  define  operator  M  for  some  large  n, 

M  =  (17) 

we  see  that 


{M[g(/0  -h  %(£,/)]}  [g{E)  -f  SgiE,t)^  dE 

- - -  (18) 

j  [giE)  +  SgiE,t)fdE 

is  in  the  form  of  equation  (12)  for  operator  M  and  its  approximate 
eigenvector  g{I^  +  Sg{E,t).  The  approximate  eigenvalue  >1„  is  accurate  to 
second  order  by  the  variational  principle  cited.  However  we  can  make 
Sg{E,t)  arbitrarily  small  by  increasing  n.  Using 

M[g(£l-f-<5g(£, /)]=;(£,/)  (19) 

and  allowing  «  oo  so  that  Sg(E,l)  -*  0,  we  rewrite  equation  (18)  as 

fy(E,/)g(£)^E 

^oo  =  —r - •  (20) 

jg\E)dE 

Comparing  this  with  equation  (13)  we  sec  that 

«(0  =  /ioc  (21) 

and  that  equation  (13)  is  variational  also. 

The  remainder 

r{E,t)  =J[E,t)  -  a{t)g{E)  (22) 

is  orthogonal  to  a(/)g(/i),  as  can  be  shown  by  multiplying  equation  (22)  by 
g(E)  and  integrating  over  energy.  If  r(E,t)  is  small  compared  to  J{E,t)  (that 
is,  if  g{E)  is  “close”  to  J[E,t)  )  then  we  expect  that  equation  (9)  will  provide 

better  estimates  of  k{t)  than  equation  (5),  given  some  uncertainty  in  the 
calculation  of  J{E,i)  due  to  error  in  electron  collision  cross-section  data.  This 
expectation  is  based  on  the  variational  property  of  equation  (13).  Uven  if 
r{E,i)  is  not  small,  it  may  be  “close”  to  another  equilibrium  energy  distrib¬ 
ution  function  g'{E)  that  is  nearly  orthonormal  to  g(£): 

jg'(£)rf£^  jg{E}g'{E)dE.  (23) 
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For  instance,  is  concentrated  at  low  energy  and  g'{E)  is  concentrated 

at  high  energy,  then  they  will  be  nearly  orthonormal.  If  equation  (23)  holds, 
the  process  may  be  repeated  to  obtain 

f  r{E,t)gV:)dE 

=  -  (24) 

J  lg'{E)fdE 

which  has  similar,  though  possibly  weaker,  variational  properties.  Clearly 
the  process  may  be  continued  in  a  manner  similar  to  a  Gram-Schmidt 
orthogonali/.ation  procedure.  The  successive  projection  cocfTicients  will 
possess  (nearly)  variational  properties  so  long  as  the  successive  g,g\ ...  are 
(nearly)  orthonormal. 

2  Solution  Method  for  Energy  Distribution  Function 


2.1  Energy  Diffusion  Equation 

To  demonstrate  the  projection  technique  described,  we  need  a  method  of 
calculating  equilibrium  and  timc-dcpcndcnt  energy  distribution  .spectra  for 
free  electrons  in  air.  Historically  most  transport  models  of  free  electrons  in 
air  have  been  based  on  a  distribution  function  which  describes  electron 
number  density  in  phase  space,  that  is,  as  a  function  of  time,  position,  and 
velocity  (or  momentum)  in  three  dimensions.  Usually  the  Maxwcll- 
Boltzmann  equation  or  some  variant  is  used  to  define  that  distribution 
function'  [chapter  2].  For  the  present  effort  wc  choose  a  simpler  formu¬ 
lation  based  on  a  distribution  function  which  depends  only  on  time  and 
electron  energy.  Derived  from  the  intcgro-differential  form  of  the  Boltzmann 
transport  equation,  this  formulation  will  be  adequate  to  develop  a  time- 
dependent  nonequilibrium  energy  distribution  of  electron  number  density 
and  to  demonstrate  the  projection  of  the  time-dependent  distribution  onto 
a  (more  or  less)  complete  set  of  equilibrium  distributions.  I'hc  projection 
will  permit  the  calculation  of  the  nonequilibrium  mobility  from  the  equilib¬ 
rium  basis. 

Because  the  position  and  direction  of  any  free  electron  in  a  uniform  gas  arc 
randomized  (except  for  drift)  after  a  few  collisions,  electron  energy  remains 
the  principal  feature  of  an  electron  swarm  (together  with  density,  which  wc 
assume  is  constant).  Certain  macroscopic  parameters  (mobility,  diffusion 
coefficient,  etc.)  can  be  calculated  from  the  electron  energy  distribution 
function.  When  the  energy  distribution  function  equilibrates  or  arrives  at  a 
steady  state,  the  equilibrium  energy  distribution  can  be  used  to  calculate 
equilibrium  macroscopic  parameters.^ 

We  will  develop  an  equation  for  the  electron  energy  distribution  function, 
the  “energy  diffusion  equation,”  with  time  and  electron  energy  as  independ¬ 
ent  parameters.  Its  terms  arc  obtained  by  balancing  energy  gains  and  losses 


^  I'rosl,  I..  S,,  and  A.  V.  Phelps,  Rotational  excitation  and  momentum  transfer  cross -sections 
for  electrons  in  lb  and  Nj  from  transport  coefficients,  Phys.  Rev.  127  (1962),  1621. 


*  Phelps,  I,.  V.,  and  I,.  C.  Pitchford,  Anisotropic  scattering  of  electrons  by  Nj  and  its  effect 
on  electron  transport,  Phys.  Rev.  A  31  (1985),  2932. 
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through  various  collisions  and  through  interaction  with  the  electric  field 
(/.cro  magnetic  field  is  assumed).  Spatial  uniformity  is  assumed,  so  spatial 
derivatives  are  omitted.  For  ambient  electric  fields  vary'ing,  say,  on  a  .scale 
of  a  tenth  of  a  nanosecond,  spatial  derivatives  can  be  neglected  if  electron 
mean  free  paths  are  much  less  than  a  few  centimeters.  In  fact  the  mean  free 
paths  are  of  the  order  of  microns  for  sea-level  air.  fhe  integro-differential 
form  of  the  Boltzmann  transport  equation  is  the  starting  point  in  the  deri¬ 
vation.  Electron  ballistic  equations  in  the  electric  field  arc  derived  and  in¬ 
corporated  into  the  transport  equation.  Reduction  of  the  transport  integrals 
to  first-  and  second-order  differential  forms  will  be  shown  for  the  case  of 
“scattering”  by  the  electric  field.  Finally  integral  terms  will  be  added  for 
scattering  by  collision  with  molecules.  Carron^  [p  9]  has  exhibited  the  en¬ 
ergy  diffusion  equation  as  obtained  from  the  equation  for  the  lowest  har¬ 
monic  of  the  electron  velocity  distribution  function.  Our  derivation  is  based 
on  a  faylor  scries  expansion  of  the  collision  integral  and  suggests  how  higher 
order  effects  might  be  estimated.  However,  our  emphasis  here  is  on  the 
projection  method  rather  than  on  the  energy  diffusion  equation  itself. 

In  the  derivation  it  will  be  assumed  that  scattering  is  isotropic,  and  therefore 
that  the  direction  cosine  of  the  emergent  electron  with  respect  to  the  ambient 
electric  field  is  uniformly  distributed  on  the  interval  (  —  I,  -t-1).  fhis  is  a  very 
good  approximation  for  clastic  scattering*  but  less  accurate  for  inelastic 
scatters.  Should  it  be  necessary,  the  probability  distribution  function  for 
angular  distribution  can  easily  be  expanded  to  include  an  additional 
Fcgendrc  polynomial  (or  more),  becoming  1  -t-  f>iX  -I-  •  ,  if  sufficient  data  for 
the  A,  are  available.  Obtaining  the  electron  drift  velocity  from  the  electron 
density  function  in  phase  space  depends  critically  upon  using  a  two-term  (or 
more)  harmonic  expansion  of  the  angular  distribution  of  the  electron  veloc¬ 
ity  (because  drift  velocity  is  closely  allied  with  the  second  term).  But  ob¬ 
taining  the  number  density  in  time-energy  space,  as  done  here,  docs  not 
depend  on  the  second  harmonic  term  for  drift  velocity  to  be  exhibited.  Both 
drift  and  diffusion  are  present  with  only  one  harmonic  term,  because  the 
harmonic  expansion  is  for  the  direction  of  the  electron  as  it  emerges  from  a 
collision,  not  for  the  direction  of  motion  of  the  electron  at  any  later  time. 

Consider  an  electron  in  air  with  electric  field  c.  The  electron's  velocity  is  u 
and  energy  L  =  \mu^.  Resolve  the  velocity  into  a  component  parallel  to  c 

M„  =  c(m  .£) 

and  a  component  normal  to  t 

\u  X  £  I 

where  c  is  the  unit  vector  parallel  to  c.  We  define  the  direction  cosine  n,  for 
«  relative  to  c  and  obtain 


»ll  =  1  I  =  wMe  (25a) 

«i  =  l«a.l  =  "\/l  •  (256) 

The  electron  will  be  accelerated  in  the  electric  field  until  it  collides  with  a 
molecule.  Assume  that  the  electric  field  is  essentially  constant  in  the  time 
between  collisions.  This  is  certainly  true  in  the  equilibrium  case  in  which  the 
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^Carron,  N.  ].,  On  the  Calculation  of  the  Electron  Energy  Spectrum  in  a  Weakly  Ionized 
Gas,  Mission  Research  Corporation  Report  MRC-R-IOSS  (30  January  1987). 


electric  field  must  be  constant  for  all  time.  Generally  this  assumption  re¬ 
quires  that  the  electric  field  change  slouiy  in  a  collision  time  (approximately 
the  reciprocal  of  the  momentum  exchange  collision  frequency).  After  a  time 
i,  u-  is  unchanged  and  u  becomes 


"ll  =  ~  (26) 

if  no  collision  has  occurred.  If  a  collision  occurs,  the  electron's  direction  is 
effectively  randomized,  on  the  assumption  of  isotropic  scattering.  fhe 
probability  of  no  collision  occurring  is  e  where  v  is  the  electron  collision 
frequency,  which  may  be  a  function  of  energy. 

As  discussed  above,  the  direction  cosine  n,  is  distributed  uniformly  on  the 
interval  (  —  I,  -f  1).  fhe  normalized  collision  probability  (unitless)  is 

dp  =  ^ve~''‘dt  dp^  (27) 


where  the  normalization  is  fdp=  1. 

In  the  usual  expression  of  the  transport  equation,  the  differential  terms  ex¬ 
press  continuity  based  on  continuous  flow  processes,  with  a  time  derivative 
plus  appropriaie  divergences  of  the  distribution  function,  fhen  all  noncon- 
tinuous  scattering  processes  are  modeled  using  collision  integrals.  However, 
continuous  or  quasi-continuous  processes  can  also  be  modeled  with  integrals 
instead  of  divergences,  and  it  is  largely  a  matter  of  convenience  which  way 
to  handle  a  process.  In  our  case,  we  use  integrals  to  model  all  scattering 
processes  and  derive  from  the  integrals  any  differential  terms  we  desire,  as 
for  scattering  by  the  electric  field,  for  example.  Thus  the  rate  of  change  of 
the  electron  energy  distribution  function /(£)  arising  from  acceleration  in  the 
electric  field  may  be  written 


MU 

di 


Coo  Too 

v{E'}/[E')dp{E',  E)  ~  v{E)AE)  dp[E,  E)  (28) 

•'£'=0 


where  dp{E',  E)  is  the  self-adjoint  differential  probability  of  an  electron  ac¬ 
celerating  (decelerating)  from  energy  E'  to  energy  E,  and  vice  versa.  The  first 
integral  accounts  for  electrons  arriving  at  energy  E  (“in-scatters"),  and  the 
second  integral  accounts  for  electrons  leaving  from  energy  E 
(“out-scatters”).  The  self-adjoint  property  of  dp(E' ,  E)  docs  not  depend  ex¬ 
plicitly  on  the  electric  field  but  docs  require  the  integral  / vdt  to  be  invariant 
under  an  exchange  of  E  and  E' .  A  sufficient  condition  for  Jvdl  to  be  invari¬ 
ant  in  this  way  is  for  the  collision  frequency  v  to  be  independent  of  energy 
over  the  range  of  energy  gain  or  loss  experienced  by  the  electron  between 
collisions.  However  when  the  electric  field  is  large,  the  energy  gain  between 
collisions  can  be  considerable  compared  to  the  initial,  unaccelerated  electron 
energy.  This  may  cause  significant  changes  in  the  collision  frequency  (or 
collision  cross-section),  but  docs  not  destroy  the  invariance  of  fvdi.  In  the 
derivation  that  lollows,  we  will  assume  that  the  self-adjoint  property  applies. 


Equation  (28)  is  exact  for  scattering  by  the  electric  field  if  the  integrals  arc 
evaluated  without  approximation.  In  general  this  is  not  possible  to  do  in 
closed  form.  One  approximation  which  allows  the  most  general  represen¬ 
tation  of  the  integrand  is  to  evaluate  the  integrals  by  Monte  Carlo;  this 
method  yields  only  numerical  results  and  entails  statistical  error  as  well,  but 
is  widely  used.  Another  approximation  which  permits  analytical  represen¬ 
tation  of  results  is  to  expand  the  integrals  as  a  Taylor  scries  in  the  energy,  a 
method  we  will  pursue.  We  want  to  expand  the  integrand  in  a  Taylor  scries 
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about  /:,  retaining  terms  only  to  second  order  in  E'  —E,  which  requires  that 
v/be  slowly  vary  ing  compared  to  dp.  However,  the  energy  dependence  of  v 
and  dp  (which  also  contains  v  )  has  not  been  specified.  Assume  for  the 
present  that  v  is  independent  of  energy.  Then  it  suffices  to  expand  /  only. 
This  causes  difficulty  when  energy-dependent  coefficients  occur  with  the  final 
differential  forms  to  be  derived,  but  a  way  out  is  ofTered  to  overcome  the 
difficulty.  Assuming  that  /  varies  slowly  compared  to  dp,  we  thus  expand 
J{E')  in  a  Taylor  series  about  E,  retaining  terms  only  to  second  order: 

AE')  ^Al^l  +  c,{E'  -E')  +  c,{E'  ~  E)\  (29) 

fhen  the  in-scatter  integral  in  equation  (28)  breaks  into  three  terms: 

Too  Too 

AE')dp{E',i-:)  =  Ai-:)\ 

Too  poo 

+  c,  {E’-E)  dp{E',E)  -t-  C2  {FJ  -  /o'  dp{E',E).  (30) 

‘'t'=T  **£'=0 

Since  dp{FJ ,  E)  =  dp{E,E'),  the  first  term  cancels  the  second  integral  in 
equation  (28),  leaving  only  the  second  and  third  terms: 

-^  =  e,v  (E'-/0dp(E',E)  +  c,v  {E' -  i:)^  dp{E' ,E).  (31) 
•’£■'=0  •'£'=0 

for  brevity  define 


=  {E' -  l^UpiE' ,i:).  (32) 

•'£'=0 

We  now  develop  a  quadrature  rule  to  approximate  the  right  hand  side  of 
equation  (31)  to  second  order,  so  that 

c^AEa)  +  bA^b)  +  cAE,)  =  c,y,  -f  (33) 

The  coefficient  a  is  a  rate  coefficient  for  scattering  electrons  from  energy  E, 
to  energy  /:*.  The  coefficient  c  is  a  rate  coefficient  for  scattering  electrons 
from  energy  F%  to  energy  /:*.  The  coefficient  />  is  a  rate  coefficient  for 
electrons  “scattering”  from  energy  to  the  same  energy.  Substituting  the 
second-order  expansion  fory(/:')  into  the  latter  equation  and  setting  Ei,=  E 
yields 

«/(/0  +.  ac^{E^  -  /:!  +  acjiEg  -  A)'  +hAi?i 
4-  cAF?i  +  CC| (/:'(.  —  /O  +  cc2{E^  —  A")'  =  C|y,  c^J-i 
Becau.se  C|  and  Ci  are  arbitrary,  this  equation  requires 

{a  +  h  +  c]Al'^  =  0 

-  ^0  +  c{f^c  -  ^]  = 

C2[a(A,-/0'  +  c(/v-/:f]  =  C2y2 
or  letting  A’„  —  A'  =  —  A  A  and  A',,  —  /:  =  A/:  and  assuming /(A")  ^  0, 
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J,-JAE 


a  = 


b  =  - 


c  = 


(A/T)^ 

jj  +  y.A/f 


2(A/:r 

Thus  the  quadrature  form  of  equation  (31), 


ol 


^  aAE  -  A/0  +  />y(/0  +  cJ{E  +  A/0. 


(34a) 

(346) 

(34c) 

(35) 


hccomcs  upon  substituting  equations  (34)  and  rearranging  terms, 

^y(/o  Ji 


at 


iy(/^+A/o-2y(/o+yi/f-A/o] 

2(A/0^ 

+  ^U[E^M-:)-AE-Am. 


Taking  the  limit  as  AE— >0  yields  the  differential  equation 

df 


1.-LJ 
dt  " 


-2  +  A  . 


dh 


(36) 


(37) 


exhibiting  terms  for  dilfusion  and  heating  of  the  energy  distribution  funetion. 
When  solutions  arc  sought  through  finitc-difrcrcncc  methods,  equation  (36) 
is  sufficient  as  it  stands,  although  differential  equation  (37)  more  succinctly 
states  the  physics  involved. 


If  equation  (29)  is  expanded  further  to  higher  powers  of  E'  —  E,  then  more 
values  of  energy  will  be  needed  in  the  quadrature  rule  expressed  by  equation 
(33).  Ihis  will  cause  the  time  derivative  ofy(/0  to  be  expanded  to  higher 
order  differences  in  (36)  and  higher  order  derivatives  in  (37).  In  fact,  by 
continuing  the  l  aylor  series  expansion  of  the  “in-scatter”  integral  it  can 
easily  be  shown  that 


df  _  y  ■/,  d‘f 

it  '■  sEt ' 


(38) 


The  importance  of  such  higher  order  terms  depends  on  their  coefficients, 
which  contain  I  he  second-order  expansion  is  central  to  the  difl'usion 
approximation  to  the  transport  equation.  Carron^  [p  12]  explains  on  phys¬ 
ical  grounds  why  the  expansion  can  reasonably  be  terminated  at  second  or¬ 
der  to  obtain  the  mobility  and  dilFu.sion  coefficient  transport  parameters. 


1; valuation  of./,  and  J2  is  tedious  and  details  arc  not  presented  here.  Using 
equations  (25)  and  (26)  we  find  that 


2  2 
9  £ 

3mv 


and 


Aq^c^E 

3mv 


-  4  4 

2q  £ 

"c  ~ 

5  m  V 


(39) 


(40) 
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where  we  have  made  the  substitution  E  =  \mu^. 

If  we  assume  that  the  collision  frequency  v  equals  the  momentum  exchange 
collision  frequency  which  is  exactly  true  for  isotropic  elastic  scattering, 
we  can  use  the  well-known  relation 


to  write 


^  = 


mv„ 


1  1  2 
./,  =  y<7^£ 

•4  =  [|  -  I  (p£/w)^]. 


(41) 


(42a) 

{A2h) 


In  general,  however,  the  (total)  collision  frequency  is  not  equal  to  the  mo¬ 
mentum  exchange  collision  frequency  when  (1)  clastic  scattering  is  no  longer 
entirely  isotropic  or  (2)  inelastic  processes  occur,  fhus,  when  dealing  with 
higher  energies  we  must  be  careful  in  our  choice  of  data  for  collision  fre¬ 
quency.  I'or  scattering  by  the  electric  field  in  nitrogen,  we  have  used  a 
combination  of  clastic  momentum  exchange  data  and  total  (or  “effective”) 
momentum  exchange  data  that  yields  energy  distributions  in  reasonable 
agreement  with  published  results,  for  oxygen  we  simply  use  total  momen¬ 
tum  exchange  data  because  of  the  paucity  of  separate  elastic  data,  fhe  main 
issue  here  is  what  collision  frequency  best  reprc.sents  the  randomization  of 
the  direction  of  the  colliding  electron,  which  is  a  complicated  topic  that  we 
shall  not  discuss  in  detail. 


■fhe  quantities  7,  and  Ji  turn  out  to  depend  on  v,  which  can  depend  on  en¬ 
ergy.  The  derivation  requires  modification  to  correctly  t.eat  this  energy  de¬ 
pendence.  A  close  inspection  of  the  integrals  shows  that  the  quantity  that 
should  be  expanded  in  a  Taylor  series  is  fjv  instead  of f.  When  fjv  is  thus 
expanded,  the  integrands  are  valid  for  energy-dependent  v.  The  elfcct  is  to 
replace  each  derivative  of / with  a  derivative  of  fjv.  first  and  second  energy 
derivatives  of  v  appear.  Because  v  varies  slowly  as  a  function  of  E,  it  is 
convenient  to  discard  .second  energy  derivatives  of  v  and  keep  only  first  en¬ 
ergy  derivatives  of  v.  This  result  is  used  by  Carron^  [p  1 3].  1  hus  the  heating 
term  becomes 


^  J-( L\ 

3m  dE  \  V  J 

and  the  diffusion  term 

d  (  E  \  ^  f  ‘  \ 

3m  BE  V  ''  dE  )  10m  BE  \  ,.2  p/.-  j 


(43«) 


(43/)) 


The  heating  term  describes  the  convection  (through  energy  space)  of 
electrons  from  lower  to  higher  energies  as  they  drift  through  the  ambient 
electric  field.  Upon  examining  the  diffusion  term,  we  observe  that  the  second 
term  in  brackets  is  smaller  than  the  first  term  by  a  factor  roughly  the  square 
of  the  electron  drift  velocity  divided  by  the  electron  speed,  fhe  second  term 
in  brackets  approximates  the  influence  of  terms  that  are  of  higher  order  than 
the  first  two  terms  in  the  usual  Legendre  expansion  of  the  velocity  distrib¬ 
ution  function.  In  derivations  of  the  velocity  distribution  function,  higher 
order  terms  than  the  second  arc  cast  away  on  the  assumption  that  the  drift 
velocity  is  much  smaller  than  the  electron  average  speed'  [p  .58  ].  fhe  as- 


sumption  that  in-flight  energy  changes  are  much  smaller  than  the  initial 
electron  energy,  upon  which  the  continuous  scattering  approximation  (cq 
(37))  is  based,  is  violated  when  the  drift  velocity  approaches  the  initial 
electron  speed.  This  occurs  in  air  at  small  electron  energies  because  of  the 
small  collision  cross-section  at  those  energies.  Tor  this  reason  we  arc  forced 
to  omit  the  second  term  in  brackets  when  dealing  with  non-idea!i/cd  air. 
Proper  inclusion  of  the  second  term  in  brackets  would  be  in  the  form  of  an 
upscattcr  probability  from  energy  level  /  to  energy  level  / -t- where  j>\. 
We  omit  this  refinement  as  outside  the  scope  of  our  investigation,  although 
possibly  deserving  of  further  examination  elsewhere.  We  retain  the  second 
term  in  brackets  in  succeeding  equations  but  omit  it  in  calculations  for  air. 


It  is  expedient  to  bring  the  coefficient  il  into  the  innermost  energy  derivative 
which  then  operates  on  the  product  Ef.  Doing  so  creates  an  additional  first 
energy  derivative  of /  which  can  then  be  combined  with  the  heating  term. 
C'ombining  the  heating  and  dilTusion  terms  in  this  way  yields  the  result 


T  2  2 

2?  £  a 

3m  dE. 


1 


d 

dE 


umtn  - 


,22 
iq  £ 

10m[v(T)]^ 


dE  v(/0 


(44) 


Using  the  same  approach  that  was  used  for  .scattering  by  the  electric  field, 
we  can  also  write  in-scatter  and  out-scatter  transport  terms  for  collision 
scattering  of  electrons  due  to  clastic  and  inelastic  cross-sections.  We  retain 
the  collision  integral  formulation  instead  of  derived  diffusion  and  heating 
differential  terms.  Inelastic  scattering  is  not  reducible  to  these  differential 
forms  because  of  large  energy  losses  that  may  occur  in  a  single  inelastic 
scatter,  preventing  the  limit  AT'  -» 0  from  being  taken  without  approxi¬ 
mation.  Although  the  difl'crcntial  approximation  is  valid  for  clastic  scatter¬ 
ing  because  the  average  energy  gain  or  loss  per  collision  is  much  smaller  than 
AT,  it  is  convenient  to  treat  both  clastic  and  inelastic  scattering  in  the  same 
way. 


The  complete  energy  diffusion  equation,  including  a  term  for  gain  of 
electrons  (from  avalanching,  ion  pair  creation,  etc.)  and  loss  of  electrons  (to 
attachment),  is 


df  2qV  e 
dt  3m  dE 


d 

dE 


,22 
3<7  c 

10m[v(/:)]^ 


d  r  y(/^d) 

dE  v(/0 


(45) 


+  p 


/*  E*—oo 

J 


u{E’)J{E)a{E,i:)dE 


•'  -  pa(/:y[ 


E*=oo 

/:*=0 


o{E,E')dE', 


where 


5(/,/:')  =  sources  and  sinks  for  electrons  of  energy  E  at  time  t, 

o{E' ,IC)  =  cross-section  to  scatter  electron  with  incident  energy  E'  into 
scattered  energy  E, 

a{E,E')  =  cross-section  to  scatter  electron  with  incident  energy  E  into 
scattered  energy  E' . 


II 


The  second  term  on  the  left  hand  side  of  equation  (45)  describes  diffusion 
due  to  collisions  in  a  uniform  electric  field.  The  third  term  describes  the 
convection  of  electrons  from  a  lower  energy  to  a  higher  energy  due  to  the 
heating  effects  of  drifting  through  a  uniform  electric  field,  fhe  right  hand 
side  describes  electron  sources  and  sinks,  and  collision  processes  between  the 
electrons  and  the  neutral  gas  molecules.  The  first  integral  accounts  for 
electrons  scattered  from  other  energies  E'  into  energy  E  (“in-scatters”),  fhe 
second  integral  accounts  for  electrons  scattered  from  incident  energy  E  to 
other  energies  E'  (“out-scatters”).  The  term  S{t,E)  may  include  processes 
proportional  to  /  and  processes  independent  of /. 

In  this  investigation  it  has  been  assumed  that  the  ioni/ation  density  of  the 
air  plasma  is  low  enough  that  only  collisions  between  electrons  and  neutral 
molecules  need  be  considered.  Thus,  electron-electron  and  cicctron-ion  col¬ 
lisions  arc  excluded  from  the  simple  model.  We  also  neglect  small-angle 
('oulomb  scattering. 

Collisions  that  must  be  included  are  elastic  collisions  and  several  kinds  of 
inelastic  collisions  such  as  excitation  of  rotational  modes,  vibrational  modes, 
and  excited  states  of  the  target  molecule,  and  ionization  of  the  molecule 
(ejection  of  an  orbital  electron).  When  used,  air  composition  will  be  taken 
to  be  oxygen  and  nitrogen  in  the  usual  proportions,  although  other  constit¬ 
uents  (water  vapor,  argon)  may  be  considered  later. 

2.2  Finite  Mesh  Model  of  the  Enerfiy  Diffusion  Equation 

We  can  solve  intcgro-diffcrential  equation  (45)  for  energy  diffusion  by  de¬ 
fining  a  finite  energy  mesh  /:„  /  =  !,  2, ...  ,  M  over  a  suitable  domain  from 
zero  to  some  large  energy  that  includes  all  electrons  of  interest.  (The  shorter 
term  finite  mesh  will  be  used  for  the  finite,  bounded  mesh  of  energies  so  de¬ 
fined.)  'l  ime  and  energy  derivatives  arc  approximated  by  finite  differences, 
and  integrals  arc  approximated  by  finite  sums.  There  are  several  finite- 
difference  approximations  to  the  differential  part  of  the  energy  diffusion 
equation,  with  varying  stability  and  accuracy  properties.  We  select  the 
Dufort-l'rankcl  scheme*  because  it  is  inherently  stable  and  simple  to  imple¬ 
ment  in  matrix  operations  on  a  computer.  Using  this  scheme,  the  differential 
terms  of  the  energy  diffusion  equation  become 


‘(Carnahan.  B.,  M.  A.  I.ulhcr.  and  J.  ().  Wilkes.  Applied  Numerical  Methods.  John  Wiley  .lad 
•Sons,  New  Y'ork  (1969),  p  451. 
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a/  2q\^  a 

dt  3m  dE 
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dE 
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?  e  a 

"*  az: 


10m[v(£)] 


3  a£' 


7AE,‘} 


AE,1) 

v(£) 
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2^  e 
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AuE  +  A£) 


v(/0 

2y(/,£) 


5m^(A/0^ 


+  ■ 


v(£- A£/2) 

At,E-^i^  ] 


[v(£  +  A£/2)]^  [v(£)]'  [v(£  -  A£/2)] 


+ 


2mA£  [  v(£ 


+  A/i)  ME-M-) 


+  A£’)  v(£-A£) 


(46) 


The  complete  finite  mesh  model  of  the  energy  diffusion  equation  is 


-,2  2 

2q  c 

2m(AE) 
+ 


c/(i  +  At,£)  -J[i  -  Af,£)] 


f  [£+ A/^f.£ 

j  v(£+A£j 


+  A/-) 


<7  e 


72) 

/(/,£  +  A£) 


2£y(/,£)  ^  [£  -  A£l/(/,£  -  AE) 


v{E)  '  v(£-A£72) 

2J[t,E)  ^  Ai,E-AE)  ] 


5m^(A£f  [  [v(£+A£/2)f  [v(£)]^  [v(£-A£/2)]^ 

£+A£)  At,E~ 
v(£ 


2  2 
«?  g 

2mA£ 

M 


f  y(/.£+A£)  At,E~AE)  ) 
I  v(£+A£)  v(£-A£')  J 


=  S(/,£)  +  pYj  tijAt,Ej)a{Ej,E^  -  p«,y(/.£i)  ^  <7(£, •,£;).  (47) 


7=1 


j 

7=1 


One  may  carry  the  sums  over  the  range  from  j  =  \  \.o  M  without  exception 
at  j  =  i,  because  the  j  =  i  terms  from  each  sum  cancel  each  other.  But  the 
notation  used  above  is  clearer. 


The  last  equation  can  be  written  in  matrix  notation  as 

M 

fit  +  Ai)  -  A/)  +  2  AyJjii)  +  So(£,-./) 


(48) 


7=1 


where 

/(/)=y(/,/a 

A,,  —  sum  of  homogeneous  terms  for  collision,  heating,  diffusion,  and 
sources/sinks,  and  where  j  and  i  are  the  column  and  row  indices,  re¬ 
spectively, 

5o(£„r)  =  inhomogeneous  .sources  and  sinks  independent  of/(such  as  ion 
pair  creation  by  radiation). 


The  matrix  elements  A,,  must  fulfill  certain  requirements  such  as  conserva¬ 
tion  of  mass,  boundary  conditions,  and  numerical  stability.  To  conserve 
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mass,  it  is  necessary  that  all  electrons  leaving  an  energy  level  arrive  at  other 
energy  levels,  or 

M 

E  4,  =  0  (49) 

(=1 

for  every  column  /,  assuming  equally  spaced  energy  levels.  The  difference 
terms  in  equation  (47)  arc  not  necessarily  mass-conservative  as  presented, 
but  calculating  A„  as  the  negative  sum  of  other  off-diagonal  column  entries 
assures  a  mass-conservative  system.  Certain  minor  changes  to  the  off- 
diagonal  entries  are  required  to  maintain  the  proper  solution  when  this  ap¬ 
proach  is  taken.  Also,  electron  attachment  cross-sections  must  be 
subtracted  from  the  main  diagonal  term.  For  numerical  stability  (and  phys¬ 
ical  rcali/ability),  it  is  necessary  that 

-1<4;<0  (50) 

or  else  more  electrons  are  removed  from  energy  level  /:,  than  exist  there,  and 

0<  Aji<  1,  ji=i  (51) 

or  else  nonphysical  “negative”  electron  densities  are  scattered  to  other  en¬ 
ergy  levels. 

Because  of  an  idiosyncrasy  of  the  Dufort-Frankcl  algorithm,  differential  er¬ 
ror  between  the  initial  conditions  J{t  —  Ai),J{t)  tends  to  be  preserved 
throughout  later  time-steps,  giving  the  appearance  of  two  solutions  differing 
by  a  constant  amount.  One  solution  occurs  at  odd  time-steps,  and  the  other 
solution  occurs  at  even  time-steps.  Replacing  y(/)  by 

after  compulation  of  each  new  y(/  -f-  At)  conveniently  unifies  the  solution. 

Elastic  Scattering 

To  solve  the  finite  mesh  equation  described  above,  the  convection  term  re¬ 
quires  the  electron  mobility  /x  as  a  function  of  electron  energy,  or  alterna¬ 
tively  the  collision  volume  as  a  function  of  electron  energy.  Flastic 
scattering  is  the  most  frequent  result  of  collision  between  conduction 
electrons  and  gas  molecules,  under  ordinary  circumstances.  The  dominant 
effect  is  to  randomize  the  direction  of  individual  electrons  (not  the  net  drift 
velocity  due  to  ambient  electric  fields)  and,  on  average,  to  transfer  momen¬ 
tum  from  the  electrons  to  the  gas  molecules.  It  is  also  possible  for  electrons 
to  gain  or  lose  small  amounts  of  energy  from  the  random  thermal  motion 
of  the  gas  molecules.  We  assume  for  our  model  of  elastic  scattering  that  the 
gas  is  monatomic  and  has  an  isotropic  Maxwellian  distribution  of  velocities, 
a  model  adequate  for  our  purposes.  The  effective  elastic  collision  cross- 
section  for  electrons  in  such  a  gas  is  equal  to  the  clastic  collision  cross- 
section  for  a  gas  without  molecular  motion,  to  a  high  degree  of 
approximation,  when  the  electrons  have  an  energy  >0.01  T  where  /'is  the 
gas  temperature^  [p  73]. 


■'Carter,  L.  L.,  and  E.  D.  Cashwcll,  Particle-Transport  Simulation  with  the  Monte  Carlo 
Method,  LSAI:RDA  Publication  TID-26607,  USAERDA  Technical  Information  Center, 
Oak  Ridge,  TN  (1975),  p  73. 


We  define  for  each  incident  electron  energy  £,  an  effective  scattering  cross- 
section  a{E„E,)  which  is  nonzero  for  j  =  i  —  \,  i,  i  +  I  and  is  zero  for  all  other 
j.  Following  a  similar  approach  to  that  used  for  scattering  by  the  electric 
field,  we  impose  constraints  to  conserve  the  first  and  second  energy  mo¬ 
ments.  (The  zeroth  energy  moment  does  not  need  special  treatment.)  Let 
the  incident  energy  E,  be  and  the  energies  E,  be  E„  Eh,  and  £  (for  nonzero 
a)  where  E,  —  Eh  —  AE  and  £  =  £*  +  AE.  Electrons  which  scatter  from  Eh  to 
Eh  do  not  require  any  cross-section,  so  we  set  ct(£4,£a)  =  0.  We  require  that 
the  effective  cross-section  be  self-adjoint,  which  is  clearly  true  for  elastic 
scattering.  The  constraints  may  then  be  written 

(£,  -  £,)a(£„£,)  -h  (£,  -  E,)aiEc,Et)  =  f  -  EMi:^  E,y  dt\  (52) 

*•0 

(£,-£,)'a(£„/4)  -h  (E,- E,fa{E„E,)  =  {E' -  E,fa{E' ,  E,)u' dE' .  (53) 

•*0 

Let  the  first  integral  be  A  and  the  second  integral  be  A  .  Suppose  the  total 
elastic  scattering  cross-section  (inversely  proportional  to  elastic  collision 
frequency)  is  independent  of  electron  energy.  This  is  often  the  case  for 
electron  energies  near  the  gas  thermal  energy,  where  electron  mobility  is  es¬ 
sentially  independent  of  ambient  electric  field.  Then 

/,  =  J  al{E,)  u'{E'  -  E,)MF,E)  dE’  (54) 


where 


/M{E‘,E)dI"7  =  self-adjoint  probability  of  scattering  from  energy  £  into 
energy  £, 

<To(/-)  =  total  elastic  scattering  cross-section  for  incident  electrons  having 
energy  £. 


It  can  be  shown  from  a  hard  sphere  scattering  model  that  free  electrons 
colliding  with  molecules  in  a  monatomic  gas  will  have  a  scattered  speed 
(magnitude  of  velocity) 


u'  =  u  — 


m 

M  +  m 


M 

M  +  m 


(55) 


w'here  M  and  U  are  the  molecular  mass  and  speed,  m,  u,  and  u'  arc  the 
electron  mass,  incident  speed,  and  scattered  speed,  C  is  the  direction  cosine 
of  the  molecular  motion  relative  to  the  bisector  of  the  incident  and  scattered 
electron  velocity,  and  ^  =  (1  +  and  where  is  the  cosine  of  the 

electron  scattering  angle. 


We  use  the  uniticss  expression 

UE,E)  dE=^  dU  (56) 

where  the  probability  of  the  target  molecule  having  speed  U  is’  [p  72] 

f^U)dU  =  -^P^  U^e-f^^'dU,  (57) 

>Jn 


and  where 

p  =  /Mi{2kri, 


IS 


k  =  Boltzmann's  constant. 

Omitting  details  of  the  integration,  we  obtain  after  discarding  small  terms 

I,^ual{E,){2kT-E)^,  (58) 

h^lual{E)EkT-fj.  (59) 

These  results  are  based  on  the  assumption  of  energy-independent  elastic 
scattering  cross-section.  If  the  clastic  scattering  cross-section  is  not  inde¬ 
pendent  of  electron  energy,  then  the  transport  integral  must  be  evaluated 
more  carefully.  Nitrogen  is  subject  to  the  Ramsauer-Townsend  effect  and 
has  an  anomalously  small  elastic  scattering  cross-section  at  low  energies 
which  varies  substantially,  leading  to  inaccurate  energy  distributions  through 
the  use  of  equation  (54).  The  epror  for  oxygen  is  less  pronounced  because 
of  lesser  variation  of  the  elastic  scattering  cross-section.  As  the  development 
of  a  solution  to  the  energy  diffusion  equation  is  of  secondary  importance  in 
this  work,  we  do  not  reduce  the  transport  integral  for  the  case  of  varying 
clastic  cross-section  but  merely  point  out  the  desirability  of  doing  so  for  the 
case  of  air. 

We  can  now'  solve 

(/',  -  E,)  W:t„EJ  +  (£,  -  E,)  c{E„E,)  =  /,  (60) 

{E,  -  E,f  a{E,,E„)  +  (E,  -  E,f  ^E„E,)  =  (61) 

for  ct(£4,  E^)  and  ct(£*,  £,)  using  E^=  E^-  AE  and  Ec=  Ei  +  AE. 

It  can  be  seen  from  the  results  for  /|  that  the  elastic  cross-section  tends  to 
scatter  hot  electrons  (relative  to  the  gas  temperature)  to  lower  energies  and 
cold  electrons  to  higher  energies,  as  expected,  until  the  median  electron  en¬ 
ergy  equals  twice  the  gas  average  temperature  kT  and  the  mean  electron 
energy  equals  \kT,  assuming  zero  ambient  electric  field.  More  than  two 
effective  cross-sections  for  clastic  scattering  arc  not  needed,  because  the  en¬ 
ergy  mesh  width  AE  is  much  larger  than  the  average  energy  change  per  col¬ 
lision. 

We  have  incorporated  these  effective  elastic  collision  cross-sections  into  the 
finite  energy  mesh  equation  described  above,  approximating  the  derivative 
terms 

and 

f  or  zero  ambient  electric  field,  the  electron  energy  distribution  relaxes  with 
time  toward  a  Maxwellian  distribution  of  energies,  as  expected,  whose  aver¬ 
age  temperature  closely  matches  the  gas  temperature.  Truncation  error  due 
to  using  an  energy  mesh  spacing  A£=0. 1  kT  and  only  60  energies  in  the 
energy  mesh  leads  to  small  errors  of  about  I  percent  in  the  final  electron 
temperature  and  0.5  percent  in  the  final  energy  distribution.  Reducing  the 
energy  mesh  spacing  to  A£'=0.05^r  and  using  120  energies  in  the  energy 
mesh  reduces  these  errors  by  a  factor  of  four  in  each  case,  suggesting  accu¬ 
racy  to  second  order  in  the  finite  mesh  equation.  The  smaller  choice  of  AE 
requires  a  smaller  At  for  stability. 

For  isotropic  clastic  scattering,  it  can  be  shown  that  the  total  clastic  scat¬ 
tering  cross-section  ah{E)  is  equal  to  the  momentum  exchange  cross-section 
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Q„.  Consideration  of  anisotropic  elastic  scatters  in  nitrogen  leads  to  changes 
not  greater  than  1  percent/  which  we  shall  neglect.  We  use  the  recent  tab¬ 
ulations  of  ^0  (total  clastic  cross-section)  for  Nj  by  Phelps  and  Pitchford* 
and  of  (momentum  exchange  cross-section)  for  O2  by  Phelps’  to  define 
o’o(E),  in  these  gases  and  in  air  composed  of  79-pcrccnt  N:  and  21-pcrccnt 
O2  by  volume. 


2.4  Unequally  Spaced  Energy  Mesh 


Up  to  this  point,  the  energy  mesh  E,  has  been  tacitly  assumed  to  be 
cquispaccd;  however  this  condition  can  be  relaxed.  It  is  helpful  to  gradually 
increase  the  energy  mesh  spacing  as  energy  increases,  so  that  fewer  mesh 
points  are  needed.  A  coarser  mesh  at  higher  energies  also  reduces  the  eficc- 
tivc  cross-section  where  it  is  largest,  thereby  permitting  a  larger  time-step  to 
be  used.  The  computation  of  the  effective  cross-sections  for  clastic  scatter¬ 
ing  is  only  slightly  changed  with  unequally  spaced  E,,  Et,  and  E^  .  We  also 
multiply 


<^(4.£'f)  by 


2iE,-E,) 

Ec-  Eg 


and  a{E)„Ea)  by 


2{E,-Eg) 

Ec-Eg 


to  account  for  unequal  energy  mesh  spacing. 


Central  differences  taken  at  E,  for  the  convection  term  and  the  diffusion  term 
become  unccntcrcd  when  an  unequally  spaced  energy  is  used,  lixpcricncc 
shows  that  unccntcring  these  terms  by  about  1  percent  of  AE  docs  not  in¬ 
troduce  important  error  to  the  solution.  If  necessary  the  convection  term 
can  be  represented  with  a  three-point  quadratic  expression  to  properly  center 
the  difference.  To  account  for  unequal  energy  mesh  spacing,  we  multiply  the 
convection  term  effective  cross-section  by  the  factors  used  for  elastic  scat¬ 
tering.  We  multiply  the  diffusion  term  effective  cross-section  by  the  squares 
of  the  same  factors. 


We  use  the  following  recursion  to  construct  an  unequally  spaced  energy 
mesh,  where  we  first  select  AEo  =  O.OSkT  and  k  =  1.01,  for  example: 


^  ^  (62a) 

A£,_,  (62h) 

AEi  =  kAE^_j  (62c) 

so  that 

AEi  =  k‘AEo  (63) 

£;.  =  A£o4rr- 


This  works  quite  well  at  all  energies.  At  low  energies  the  spacing  is  nearly 
AEn,  and  at  high  energies  the  spacing  is  nearly  (k— 1)E.  For  zero  electric 
field,  the  error  in  the  energy  distribution  function  at  the  high  end  of  the  en- 


*  Phelps,  A.  V.,  and  L.  C.  Pilchford,  Anisotropic  Scattering  of  Electrons  by  N2  and  its  Effects 
on  lilectron  Transport:  Tabulations  of  Cross  Section  and  Results,  JlEA  Information 
Center  Report  No.  26,  University  of  Colorado,  Boulder,  CO  (I  May  1985),  p  14. 


’  Phelps,  A.  V.,  Tabulations  of  Collision  Cross  Sections  and  Calculated  Transport  and  Re¬ 
action  Coclficicnts  for  F-lcctron  Collisions  with  O2,  JILA  Information  Center  Report  No. 
28,  University  of  Colorado,  Boulder,  CO  (1  September  1985),  p  10. 
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ergy  mesh  is  less  than  1  percent  for  the  mesh  constants  selected,  until  the 
distribution  is  reduced  in  amplitude  by  over  10  decades  at  the  high-energy 
end  of  the  mesh. 

For  inelastic  scattering,  variable  energy  mesh  spacing  can  lead  to  systematic 
error  and  therefore  must  be  used  with  caution.  We  have  tried  several  ways 
of  correcting  the  cfrcctivc  inelastic  cross-sections  to  account  for  variable  en¬ 
ergy  mesh  spacing  without  complete  success.  As  a  minimum  correction,  the 
effective  cross-sections  for  inelastic  scattering  arc  divided  by  k.  For  example, 
results  for  the  energy  distribution  arc  low  by  about  15  percent  for  energies 
above  several  electron-volts,  using  k  =  1.01,  because  of  dominance  by  rap¬ 
idly  changing  inelastic  cross-sections  at  energies  of  several  electron-volts. 
We  must  use  a  variable-spaced  mesh  when  covering  thermal  energies  to 
ionizing  energies  (above  roughly  14  cV)  because  approximately  10,000 
equispaced  energies  would  be  needed.  Dense  matrices  of  such  size  arc  too 
large  for  available  computers. 

2.5  Exact  Solutions  and  Equilibrium  Finite  Mesh  Solutions  Compared 

Exact  solutions  for  the  equilibrium  electron  energy  distribution  arc  known' 
[pp  71-75]  when  scattering  is  limited  to  clastic  scattering  and  (1)  electron 
drift  velocity  is  much  less  than  the  mean  molecular  speed  (.Maxwellian  dis¬ 
tribution)  or  (2)  electron  drift  velocity  is  much  greater  than  the  mean  mo¬ 
lecular  speed  and  the  clastic  cross-section  is  independent  of  energy 
(Druyvcstcyn  distribution).  In  the  first  case  the  time  between  collisions  (re¬ 
ciprocal  of  collision  frequency)  is  independent  of  electron  energy;  in  the  sec¬ 
ond  case  the  pathlcngth  between  collisions  (or  collision  cross-section)  is 
independent  of  energy. 

For  zero  ambient  electric  field,  the  Maxwellian  distribution  of  electron  en¬ 
ergies  is 

y(/0  =  (65) 

For  large  ambient  electric  field,  the  Druyvcstcyn  distribution  of  electron  en¬ 
ergies  is 

(66) 

where  is  independent  of  electron  energy. 

It  is  easy  to  show  that  the  energy  mesh  equation  (cq  (47))  correctly  models 
these  distributions  by  applying  the  principle  of  detailed  balance  to  the  energy 
distribution  function.  For  our  case  of  a  tridiagonal  transition  matrix  where 
transitions  only  occur  between  adjacent  energies,  detailed  balance  states  that 
in  the  equilibrium  limit  the  rate  of  electrons  going  from  energy  E  to  energy 
E  4-  SE  equals  the  rate  of  electrons  going  from  energy  E  -f  SE  to  energy  E. 
Assume  an  equispaced  energy  mesh.  Define  the  rate  of  upgoing  electrons 
as  R*{E\f{E)  and  the  rate  of  downgoing  electrons  as  R  {E  +  Af^E  +  AE). 
Then  detailed  balance  gives 

AE)  ^  R~{E+AE) 

AE  +  AE)  ^ 

Take  the  case  of  the  Maxwellian  distribution,  where  c  =  0  and  therefore 
7i  =  =  0.  From  previous  derivations,  we  have 
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(OM 

R-{E  +  - - ^—r  {[\  -  /',  A/0  (68fc) 

2{AE) 

where  the  primed  quantities  are  evaluated  at  /:  +  A/:  and  the  corresponding 
unprimed  quantities  are  evaluated  at  E.  Then 

//(/:  + A/Q  _  (T,  -/',A/0 

//(/O  (/2  +  /,a/0  ■  ^  ^ 

Assuming  ol{E)  ~  qIu{E)  where  ^  is  independent  of  energy,  and  dropping 
second  and  higher  powers  of  AE,  we  obtain  after  some  manipulation 

i^AIL  +  AIL 

R-{E+AE)IR\E)  ^  (70) 

*  ~  2kT  '^~ir 

which  equals,  to  the  order  of  AE  retained,  the  exact  solution 


Similarly,  for  high  electric  fields  such  that  kT  4:  E,  and  also  neglecting  the 
high-order  correction  diffusion  term,  the  downscattcr  and  upscatter  transi¬ 
tion  rates  between  energy  levels  E  +  AE  and  E  arc 

. . .  29V  [E+M-yiu+iii-:) 

R  (/:  -I-  AE)j{E  -f  A/0  - - 5 - A.ri~>\ 

3m(A/0^  v{t  +  AEI2) 


AE+Af-) 
2mAE  v{E+AE} 


jrf-  v(E  +  A/02) 
mpuaiAE 

/T  r'  ,  A 


2M{AIC)‘ 


■AE+AE) 


R{IW^  =  -  ""  -r^/O  (72/2) 

3ot(A/0  2mA/:  v(/0  2M{AEf 

where  terms  involving  A; 7' have  been  neglected.  We  now  apply  equation  (67). 
Using  moo  =  Qm  and  v  =  NuQ„,  after  some  nrianipulation  wc  arrive  at 

AE-AE)  \F 

J^-^=X-^^2BEAE  (73) 


3m  (  ^Qn 


where 


which  agrees  to  first  order  with  the  Druyvesteyn  distribution 


y(/0  oc  ^fE 


In  order  to  establish  convergence  and  consistency  of  our  numerical  solution 
to  the  energy  mesh  equation  (eq  (47)),  w'c  numerically  calculate  equilibrium 
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solutions  for  the  (65)  and  (66).  All  elements  Aj,  of  matrix  A  in  equation  (48) 
must  be  smaller  than  unity  in  magnitude  for  the  equation  to  converge  nu¬ 
merically.  Otherwise,  numerical  instabilities  occur.  The  elements  on  the 
main  diagonal  and  first  upper  and  lower  subdiagonals  contain  terms  pro¬ 
portional  to  the  square  of  the  electric  field  c.  These  terms  become  very  large 
as  the  electric  field  becomes  large.  The  time  stepsi/e,  a  multiplier  of  these 
terms,  must  be  reduced  to  very  small  values  to  keep  the  corresponding  ele¬ 
ments  smaller  than  unity  and  maintain  numerical  stability.  Solutions  for  the 
energy  distribution  function  may  take  several  microseconds  to  relax  toward 
equilibrium.  If  the  time  stepsize  is  reduced  to  1  fs  (1  fs  =  0.001  ps)  to 
maintain  stability,  for  example,  thousands  of  hours  of  computer  running 
time  will  be  needed  to  find  an  equilibrium  solution.  However,  the  number 
of  iterations  of  the  time-advancement  algorithm  can  be  reduced  drastically 
by  increasing  the  number  of  time-steps  the  energy  distribution  function  is 
advanced  at  each  iteration.  I'his  is  done  by  constructing  a  single  matrix 
which  advances  the  distribution  function  vector  many  time-steps  at  a  single 
multiplication. 

Omit  the  inhomogeneous  term  So{E„t)  from  equation  (48)  and  define  the 
notation  y;;  =  ^(mA/).  Then  the  equation  becomes 

fn=fn-l^  ^fn-U  (74) 

fhe  averaging  process  of/i  and  f„  2  to  reduce  differential  error  in  the  initial 
conditions becomes 

=\{fn+fn-2l  (75) 


where is  subsequently  used  in  place  of Assume  initial  conditions 
/o  =yi  =  0.  Then  f„  is  defined  recursively  by 


fn  =  (76) 

f'n-X  (77) 


Introduce  matrices  B„  and  C„  such  that  f„  =  B„fo  and i  =  C„fa.  Recursions 
for  these  matrices  can  easily  be  derived,  based  on  equations  (76)  and  (77): 
We  can  control  the  growth  of  numerical  roundoff  error  by  normalizing  the 
determinant  of  matrix  B  to  unity  at  each  recursion.  I'hus  with  the  above 
initial  conditions  and  defining  I  =  identity  matrix,  we  have  the  recursion 

B,  =  1,  C,  =  I, 

B;^,  =AB„-t-C„,  (78) 


B„ 


B' 


n+\ 


IIBW,  II  ’ 

^n+i  =  y(B„+i  +  C„). 


(79) 

(80) 


Now  that  we  have  B,  for  any  n,  where  f„  =  B„_/o,  we  can  advance  from  f„  to 
yi,  by  using  yL,  =  B„fn.  This  requires  only  the  additional  assumption  that 
=  fn,  which  is  a  good  approximation.  The  function /  will  relax  toward  the 
equilibrium  state  for  any  reasonable  choice  of  initial  conditions,  and 
equation  (80)  will  guarantee  that  differential  errors  diminish  (being  roughly 
halved  at  each  application).  Taking  n  =  10,  for  example,  reduces  differential 
error  by  about  a  factor  of  a  thousand.  Thus,  the  function  / can  be  advanced 
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I'initc  mesh  equilibrium  results  for  Maxwellian  distribution:  Solid  line  is  finite  mesh  result,  dashed 
line  is  Maxwellian  distribution  for  electron  temperature  7’ such  that  kr=  0.025  cV,  and  dolled  line 
is  magnitude  of  the  difTcrcncc.  The  two  distributions  have  been  matched  at  the  maxima  to  eliminate 
normalization  error.  Arbitrary  units  are  used  for  the  dependent  axis.  Only  elastic  scatters  arc  con¬ 
sidered. 


in  leaps  of  2"  x  10  time-steps  by  (1)  calculating  Bio  using  the  above  recursion, 
(2)  calculating  matrix 

^  (B,of  (81) 


by  using  the  recursion 

Do  =  B,o,  =  (82) 

and  (3)  multiplying /by  D„  for  each  leap.  For  the  appropriate  choice  of  m 
in  evaluating  matrix  D„,  equilibrium  solutions  can  be  calculated  rapidly  even 
for  large  values  of  the  electric  field.  Additionally,  an  elementary  result  of 
linear  algebra  shows  that,  as  m  oo  ,  the  columns  of  D„  each  approach  the 
eigenvector/,,  =  g(/:,)  multiplied  by  a  dilTerent  scalar  factor  for  each  column 
(hence  the  infrequent  usage  “eigencolumn”  for  “eigenvector”). 

To  test  the  finite  mesh  equation,  we  have  obtained  equilibrium  solutions  by 
iterating  the  time  leap  (for  «  =  10  and  m  =  25)  until  the  energy  distribution 
became  unchanging,  for  mesh  constant  A  =  I  for  the  Maxwellian  case  and 


21 


ENERGY  (aV) 


Kinitc  mesh  equilibrium  results  for  Druyvcstcyn  distribution:  Solid  line  is  finite  mesh  result,  dashed 
line  is  Druyvesieyn  distribution  for  electron  mean  energy  equal  to  100  limes  gas  thermal  energy, 
where  k  T  =■  0.025  cV,  and  dotted  line  is  magnitude  of  the  difference.  1  he  two  distributions  have  been 
matched  at  the  maxima  to  eliminate  normalization  error.  Arbitrary  units  are  used  for  the  dependent 
axis.  Only  elastic  scatters  are  considered. 


k  =  1.01  for  the  Druyvcstcyn  case,  using  =  k  l'l20  =  0.00125  cV.  I'or  the 
Maxwellian  case,  the  elastic  cross-section  cto  was  made  to  be  constant.  For 
the  Druyvcstcyn  case,  ctJ  was  made  to  vary  as  l/w{F)  except  at  the  lowest 
energies,  where  al  was  made  nearly  constant  to  avoid  the  singularity  at 
li  =  0.  Also,  the  high-order  correction  dilTusion  term  was  neglected.  .Mean 
electron  temperatures  equal  to  kT  (Maxwellian  distribution)  and  lOOAT 
(Druyvcstcyn  distribution)  were  used.  A  lower  boundary  condition  was  ap¬ 
plied  which  rcquircsy(/;,)  =J{Ii2)lj2  .  Calculated  results  for  the  Maxwellian 
distribution  and  error  compared  to  the  exact  solution  arc  shown  in 
Figure  1  on  page  21.  Similar  results  for  the  Druyvcstcyn  distribution  arc 
shown  in  I'igure  2,  As  can  be  seen,  the  energy  mesh  equation  gives  satis¬ 
factory  results  even  many  orders  of  magnitude  down  from  the  distribution 
maximum.  Xoticablc  error  occurs  at  the  lowest  energies  because  of  the  large 
spacing  of  the  energy  mesh  compared  to  the  electron  energy. 
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2.6  Inelastic  Scattering 


Inelastic  scatters  of  incident  electrons  by  gas  molecules  can  involve  large 
changes  in  energy,  in  contrast  to  small  changes  caused  by  clastic  scattering 
and  scattering  by  the  ambient  electric  field.  For  small  changes,  the  zeroth, 
first,  and  second  energy  moments  of  scattering  probability  (integrating  over 
the  entire  energy  spectrum)  can  be  represented  accurately  using  distribution 
function  values  at  only  three  adjacent  energies  in  the  energy  mesh.  We  have 
shown  how  such  a  three-point  scheme  approximates  the  sum  of  a  convective 
first  derivative  and  a  dilTusive  second  derivative  with  the  derivatives  evalu¬ 
ated  at  the  scattered  electron  energy,  which  is  equivalent  to  expanding  the 
energy  distribution  function  to  second  order  in  the  neighborhood  of  the 
scattered  electron  energy  to  estimate  the  “in-scatter”  collision  integral,  fo 
preserve  accurately  the  same  first  three  energy  moments  of  scattering  prob¬ 
ability  when  large  changes  in  energy  occur,  it  is  necessary  to  expand  the  en¬ 
ergy  distribution  function  to  second  order  in  the  neighborhood  of  the 
incident  electron  energy.  Because  the  energy  change  is  fixed  for  each  sepa¬ 
rate  cxcitation/de-excitation  process  (except  for  ionization),  the  integral  of 
the  energy  moment  of  the  process  cross-section  reduces  to  a  product.  That 
is,  the  probability  of  the  electron  scattering  from  E'  to  £  is  a  delta  function, 
5(7:' -/O-  Thus, 


f  uiE')a,il7)iI‘:  -  Ef  df  =  (SE.ruiE,  SE.Hl-'i  + 

Jn 


(83) 


where 

o*(/i)  =  Ath  process  cross-section, 

SEj,  =  kth  process  energy  loss, 

/:,  =  /th  energy  in  the  energy  mesh. 

Although  the  scattered  energy  will  always  be  a  mesh  energy  /:„  the  incident 
energy  will  generally  lie  between  two  mesh  energies  Ej  and  E^,,.  The  cross- 
section  given  for  the  incident  energy  must  therefore  be  approximated  by 
equivalent  cross-sections  at  the  two  mesh  energies  E,  and  7:^,,.  The  equiv¬ 
alent  cross-sections  arc  chosen  in  such  a  way  that,  when  combined,  they 
conserve  the  actual  cross-section  and  its  first  energy  moment.  Because  only 
two  mesh  points  arc  used,  an  error  is  introduced  in  second  and  higher  energy 
moments  which  diminishes  in  relative  importance  for  larger  energy  losses. 

For  N2  rotational  excitation  we  use  the  two-term  approximation  proposed 
by  Phelps  and  Pitchford*  (upper  table  therein).  As  explained  by  Goldstein, 
this  two-term  approximation  incorporates  a  continuous  scattering  term  to 
represent  rotational  excitation  at  electron  energies  less  than  0.8  cV  and  a 
single-level  excitation  term  with  0.02-eV  energy  loss  to  represent  rotational 
excitation  near  the  2-cV  resonance  of  the  Nj  molecule.  This  two-term  ap¬ 
proximation  was  adjusted  by  its  originators  to  reproduce  the  diffusion  coef¬ 
ficient  and  mobility  parameters  for  the  electron  distribution  function,  fhe 
continuous  approximation  used  includes  only  a  first  derivative  term,  instead 
of  first  and  second  derivative  terms  as  recommended  by  Carron  recently,*  but 
reproduces  the  principal  transport  coefficients  in  the  energy  range  where  ro¬ 
tational  excitation  is  important.  It  should  be  pointed  out  that  the  collision 


I®  Goldstein,  B.,  A  Summary  of  Rotational  and  Vibrational  Cross  Sections  in  Nj,  Mission 
Research  Corporation  Report  MRC-R-IOS?  (26  January  1987). 
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operator  described  by  Goldstein  (his  equation  (1),  drawn  from  Frost  and 
Phelps^)  and  also  the  two-term  approximation  (his  equation  (11))  are  valid 
for  the  steady-state  coupled  equations  that  I'rost  and  Phelps  obtained  by 
expanding  the  electron  distribution  function  in  a  harmonic  expansion  of 
Legendre  polynomials,  and  are  not  the  correct  collision  operators  for  the 
time-dependent  Bolt/mann  equation  as  a  function  of  velocity  or  energy.  I'or 
example,  the  coefTicients  for  inelastic  collision  terms  in  Goldstein's  ex¬ 
pression  of  the  collision  operator  must  be  divided  by  energy  loss  and  multi¬ 
plied  by  electron  speed  to  be  used  with  the  time-dependent  Boltzmann 
equation  described.  For  O2  rotational  excitation,  we  use  the  single-level  ap¬ 
proximation  proposed  by  Phelps.’ 

Vibrational  excitation  cross-sections  for  N2  arc  taken  from  Phelps,*  which 
includes  tabulations  for  transitions  from  ground  state  (/  =  0)  to  each  of  the 
first  eight  excited  states  (i>"=l  to  8).  Ground  state  was  the  only  state  as¬ 
sumed  to  be  populated  among  the  neutral  molecules.  Similarly,  O2  cross- 
section  tabulations  for  transitions  from  ground  state  to  the  first  four  excited 
vibrational  states  were  used,’  including  both  low-cncrgy  and  9-cV  resonance 
data.  Electronic  excitation  and  ionization  cross-section  tabulations  were 
drawn  from  the  same  sources. 

I'or  rotational,  vibrational,  and  electronic  excitation  processes,  the  energy 
lost  by  the  incident  electron  is  well-defined  and  fixed  for  each  transition. 
However,  for  ionization  the  residual  kinetic  energy  must  be  partitioned  be¬ 
tween  the  scattered  electron,  the  ejected  electron,  and  a  possible  excited  state 
of  the  target  molecule.  Cross  sections  for  ionization  excitation  of  N2  are 
available"  but  are  scarce  for  O2.  Some  results  are  likewise  available  on  the 
energy  distribution  of  the  ejected  electron.  For  the  present  study,  we  will 
simply  assume  a  ground  state  N2  product  and  divide  the  remaining  kinetic 
energy  equally  between  the  scattered  electron  and  the  ejected  electron. 

fwo-body  and  three-body  attachment  cross-sections  are  used  for  molecular 
oxygen.’  fhe  term  “equilibrium  energy  distribution”  in  oxygen  or  air  (or 
other  attaching  gases)  requires  special  interpretation,  because  all  free 
electrons  ultimately  attach.  We  interpret  the  equilibrium  energy  distribution 
for  an  attaching  gas  in  the  absence  of  electron  sources  to  be  the  limit  of  the 
energy  distribution  as  the  electron  density  goes  to  zero.  In  general,  such  a 
limit  will  exist.  When  electron  sources  exist,  such  as  ionization  sources  or 
avalanching  (breakdown)  sources,  different  interpretations  arc  required.  I'or 
a  constant  source  of  electrons  such  as  time-independent  ionizating  radiation, 
the  electron  density  will  tend  toward  a  finite  limit  and  the  corresponding 
equilibrium  energy  distribution  will  be  defined  as  for  nonattaching  gases. 
For  an  increasing  source  of  electrons  through  avalanching,  gas  molecules 
would  ultimately  become  stripped  of  electrons  and  electron  density  limited. 

I  lowcvcr,  we  would  be  more  interested  in  a  quasi-equilibrium  state  where 
electron  density  growth  is  exponential  and  electron  energy  has  attained  some 
temporarily  stationary  distribution,  fhis  quasi-equilibrium  distribution  is 
the  interpretation  we  would  use  in  such  a  case. 
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"  Wadzinski,  l(.  T.,  and  J.  R.  jasperse,  l.ow  I'.ncrgy  Idcctron  and  Photon  Cross  Sections  for 
O,  N2.  and  O2.  and  Related  l)ata.  Air  Force  (leophysics  Laboratory  (PIIY),  Report 
Af  Gl.  TR  82  0008  (4  January  1982). 
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Figure  J.  Finite  mesh  results  and  first  two  Legendre  cocflicients  compared;  Solid  line  is  Finite  mesh  result, 
dashed  line  is  first  Legendre  coefTicicnt  of  Phelps  and  Pitchford,  dash-dotted  line  is  second  Legendre 
cociricienl,  and  dotted  line  is  sum  of  two  Legendre  cocITicicnts.  The  Legendre  cociricients  have  been 
normalized  so  that  the  maximum  of  the  sum  is  unity.  The  maximum  of  the  finite  mesh  result  is  also 
unity.  The  finite  mesh  result  has  been  divided  by  the  square  root  of  the  energy  to  conform  to  the 
l.cgcndrc  coefficients.  All  curves  were  calculated  for  c/A'  =  100  Td  in  nitrogen. 


2.7  Equilibrium  Energy  Distribution  Compared  to  Previous  Results 

Figure  3  compares  equilibrium  energy  spectrum  results  of  the  finite  mesh 
equation  with  our  rough  digitizations  of  Phelps  and  Pitchford's  equilibrium 
energy  spectrum'^  for  pure  molecular  nitrogen.  1'he  spectra  were  calculated 
for  c//V=  100  Td  (I  Townsend  or  Td  =  1  x  10  ’’  V-cmVmolecule)  where  N 
is  molecular  density  and  c  is  electric  field  strength.  The  complete  set  of  ine¬ 
lastic  cross-sections  for  Ni  was  used  for  the  finite  mesh  equation.  The  energy 
mesh  was  constructed  for  k=  1.015  and  A/ii  =  kTjlO,  except  that  a  maxi¬ 
mum  mesh  spacing  oC2kT was  permitted.  The  finite  mesh  result  is  compared 
with  the  first  two  Legendre  coefficients  of  a  six-term  solution  for  the  electron 
energy  distribution  function.  The  sum  of  the  first  two  Legendre  coefficients 


'^Pitchford.  L.  C.,  and  A.  V'.  Phelps,  Comparative  calculations  of  electron -swarm  properties 
in  \2  at  moderate  A'/ A' values,  Phys.  Rev.  A  25  (1982),  540-554. 
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compares  well  with  the  finite  mesh  equation  result.  The  latter  correctly  re¬ 
produces  the  “bump”  at  about  2  cV  caused  by  the  second  Legendre  coefil- 
cient.  The  sharp  drop  in  the  finite  mesh  result  at  energies  near  zero  agrees 
with  a  Monte  Carlo  calculation  done  by  Phelps  and  Pitchford,  although  our 
drop  is  steeper.  Wc  attribute  the  small  number  of  electrons  near  zero  energy 
to  the  exceptionally  small  collision  cross-section  in  this  energy  regime  for 
nitrogen  (the  well-known  Ramsaucr-Townsend  eficct),  allowing  electrons  to 
accelerate  to  significantly  larger  energies  in  the  electric  field  between  colli¬ 
sions.  The  high-energy  tail  is  smaller  for  the  finite  mesh  result  than  for  the 
Phelps  and  Pitchford  distribution,  but  rough  agreement  is  observed.  The 
average  energy  of  the  finite  mesh  result  is  2.0  eV,  compared  to  2.2  cV  for  the 
Phelps  and  Pitchford  distribution.  Noticeable  error  is  seen  at  various  ener¬ 
gies,  however,  which  underscores  the  possible  importance  of  a  variational 
method  which  might  reduce  errors  in  the  energy  distribution  to  second  order 
in  the  computation  of  swarm  parameters. 


s/A  (Td) 

Mesh  (eV) 

Average  energy  factor 

Drift  velocity  (m/s)  | 

calculated 

measured 

calculated 

measured 

0 

^kT 

20^' 

0.77 

1.00 

0.00 

0.00 

0.01 

^kT 

20^  ‘ 

1.21 

— 

3.5  X  10^ 

4.0  X  10^ 

0.1 

4.3 

1.8 

1.1  X  W 

2.5  X  10’ 

0.3 

^kT 

6.4 

4.8 

2.3  X  W 

4.0  X  10’ 

0.3 

\kT 

6.5 

4.8 

2.3  X  10’ 

4.0  X  10’ 

1.0 

\kT 

13. 

12. 

4.6  X  10* 

5.0  X  10’ 

10. 

ikT 

28. 

36. 

2.7  X  10* 

2.0  X  10* 

100. 

\kT 

63. 

72. 

1.5  X  10* 

1.0  X  10* 

Table  I.  Calculated  and  Measured  Swarm  Parameters  in  Nitrogen  Compared;  Approximate  swarm  param¬ 
eters  from  measured  data  and  swarm  parameters  from  finite  mesh  equation.  Swarm  parameters  are 
averaged  over  the  electron  energy  distribution,  electron  energy  factor  is  average  electron  energy  di¬ 
vided  by  average  molecular  energy  ^kT. 


2.8  Equilibrium  Transport  Properties  Compared  to  Measured  Data 

Wc  calculated  ensemble  (swarm)  mobility,  drift  velocity,  and  average 
electron  energy  in  nitrogen  for  a  small  number  of  c//V  values  from  zero  to 
100  I’d,  using  equilibrium  solutions  of  the  finite  mesh  equation.  Results  are 
shown  in  Tabic  1.  Because  of  significant  errors  caused  by  variable  energy 
mesh  spacing  when  used  with  inelastic  cross-sections,  two  constant  mesh 
spacings  were  used.  A  fine  mesh  (A£’=  ^kT)  was  used  for  e/jV  from  zero  to 
0.3  Td,  and  a  coarse  mesh  (A/''=  \kl')  was  used  for  e//V  from  0.3  to  100  fd, 
with  an  overlap  at  0.3  Td.  Average  electron  energy  is  given  in  units  of 
\kT. 

Inaccuracy  due  to  energy  dependence  of  the  elastic  scatter  cross-section  in 
nitrogen  can  be  seen,  as  indicated  earlier  in  the  discussion  of  clastic  scatter¬ 
ing.  The  average  energy  should  tend  toward  unity  as  the  electric  field  goes 
to  zero,  but  a  different  limit  is  observed.  This  is  strictly  due  to  the  energy 
dependence  of  the  clastic  scatter  cross-section,  as  established  by  additional 
calculations  based  on  a  fixed  cross-section.  Comparisons  with  approximate 
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values  taken  from  published  measured  data'^  range  from  poor  to  good.  Ac¬ 
curacy  of  the  approximate  values  for  the  data  is  not  more  than  10  or  20 
percent,  given  the  scatter  of  data  and  the  difficulty  of  digitization  by  in¬ 
spection  of  the  original  graphs.  The  two  finite  mesh  results  at  0.3  Td  for  the 
fine  and  coarse  energy  meshes  agree  well  with  each  other.  Finite  mesh  re¬ 
sults  are  in  good  agreement  with  data  at  1,  10,  and  100  Td.  Agreement  is 
relatively  poor  at  0,  0.01,  and  0.1  Td  because  of  the  varying  elastic  cross- 
section. 

If  our  objective  were  to  calculate  accurately  the  swarm  parameters  directly 
from  the  electron  energy  distribution,  we  might  strive  for  more  accuracy  than 
these  results  indicate.  However,  the  objective  of  this  work  is  to  demonstrate 
and  evaluate  a  method  which  may  give  good  estimates  of  swarm  parameters 
'n  the  presence  of  various  errors  introduced  by  the  cross-section  set  and  by 
the  computation  of  the  electron  energy  distribution  function.  Indeed,  as  will 
be  shown,  the  method  trivially  reproduces  experimentally  determined  equi¬ 
librium  swarm  data.  The  most  important  test  of  the  method  will  be  calcu¬ 
lation  of  nonequilibrium  swarm  data,  however. 

3  Analysis  of  Equilibrium  Projection  Method 


3.!  Developing  the  Basis  of  Equilibrium  Energy  Distributions 

We  now  desire  to  construct  the  basis  0(£,  A,)  used  in  equation  (6).  For 
convenience,  the  parameter  A.  is  replaced  by  the  equivalent  mapping 
c,IN),  where  A,  is  identified  with  the  density-normalized  electric  field  c, 
in  gas  having  molecular  density  N.  For  concreteness,  we  specify  the  gas  to 
be  air  as  previously  constituted.  We  restrict  the  basis  to  finite  values  of  /, 
letting  i  run  from  1  to  n.  The  c,  are  chosen  so  that  the  corresponding  <b  are 
nondegenerate  and  complete  in,  for  example,  the  finite  space  E  =  Ej.  That 
is,  a  suitably  bounded  function  F(Ej)  can  be  approximated  as 

n 

y=  1.  Af.  m 

1=1 

This  description  arises  from  fitting  the  continuous  function  F{E)  at  the 
points  Ej  with  a  finite  sum  of  linearly  independent  basis  functions 
(each  multiplied  by  a  coefficient /»,).  If  M  >  /?,  the  resulting  linear 
system  of  equations  is  ovcrdctcrmincd  and  a  solution  is  obtained  typically 
by  minimizing  a  measure  of  the  error  of  the  approximation.  If  M  =  /i  the 
system  can  be  solved  exactly  as 

n 

F{Ej)  =  j=  1,...,  n.  (85) 

/=! 

In  our  application  the  functions  F  and  4>  represent  continuous  functions  of 
electron  energy  E.  Thus,  solving  the  latter  equation  so  that  F  is  fitted  exactly 
at  n  energies  Ej  does  not  guarantee  acceptable  behavior  of  the  fit  at  inter¬ 
vening  energies  (or  energies  outside  the  range  of  E^.  In  that  which  follows 
we  must  be  aware  that  a  poor  fit  at  these  other  energies  may  drastically  alter 
the  calculation  of  collision  volume  (for  example)  for  the  fitting  function 


DuUon,  J.,  A  survey  of  electron  swarm  data,  J.  Phys.  Chem.  R^.  Data,  4  (197S),  S77-8S6. 
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shown  on  the  right  hand  side  of  equation  (85).  Conversely,  calculating  a 
nonequilibrium  collision  volume  using  a  measurement  of  equilibrium  colli¬ 
sion  volume  implies  knowledge  of  the  distribution  at  all  energies,  an  impli¬ 
cation  not  consistent  with  the  restriction  M  =  n.  This  means  that  the  use 
of  a  measurement  of  equilibrium  collision  volume  may  yield  bad  results  if 
no  consideration  is  given  to  the  goodness  of  fit  at  intervening  energies. 

A  trial  basis  may  be  constructed  from  any  set  of  n  distil >  values  of  cJN. 
However,  because  compulations  with  finite  precision  will  be  used,  it  will  be 
helpful  to  select  values  that  promote  linear  independence  among  equilibrium 
energy  distributions  derived  therefrom.  This  can  be  done,  for  example,  by 
scattering  the  values  over  the  range  of  interest,  say,  from  zero  to  a  selected 
maximum  value,  fhe  values  may  be  unequally  spaced,  if  desired. 

The  coefficients  p,  describe  the  spectral  projection  of  /'(/:,)  onto  the  (ft  basis 
and  are  not  uniquely  defined  while  the  basis  vectors  arc  linearly  dependent. 
In  matrix  notation  with  row-column  order  of  subscripts,  so  that  matrix 
(RX,  =  f./A')  and  vectors  (p),  =  p,  and  (f),  =  /•'(£)),  we  have 

f=Rp.  (86) 

Each  column  of  matrix  R  is  a  basis  vector  (ft  for  some  choice  of  c,/.V.  Be¬ 
cause  the  collision  volume  K  for  an  ensemble  is  linearly  dependent  upon  the 
energy  distribution /(/:■)  (sec  cq  (4)),  we  can  write  the  ensemble  total  collision 
volume  in  the  spectral  form 

n 

L,al=Y,kPi  (87) 

(=1 
A 

where  K,  is  the  collision  volume  for  an  ensemble  whose  equilibrium  energy 
distribution  corresponds  to  a  choice  tJN  for  electric  field  divided  by  molec- 

A 

ular  density.  In  matrix  notation  where  (k),  =  K„  we  write  this  as 

L,ai=y^''p-  (88) 

Solving  equation  (86)  for  p  and  substituting  in  the  last  equation  gives 

(89) 

At  this  stage,  several  methods  of  obtaining  unique  projection  coclTicicnts  p, 
must  be  considered.  One  method  is  to  take  M  =  n  and  solve  equation  (85) 
for  the  projection  coefficients.  However,  this  method  is  not  truly  consistent 
with  our  desire  to  obtain  a  smooth,  close  fit  to  the  nonequilibrium  distrib¬ 
ution  at  intervening  energies.  Hence  we  desire  to  make  M  P  n  and  take  ac¬ 
count  of  as  many  energies  as  practical.  The  coclTicicnts  p,  may  take  on 
positive  or  negative  values  as  needed  to  manage  the  fit  at  all  the  E„  con¬ 
tributing  to  a  numerically  unstable  solution.  Another  method  is  to  impose 
a  constraint  that  p,  >  0  for  each  /,  and  find  an  approximate  solution  to 
equation  (86)  that  satisfies  this  and  perhaps  other  constraints.  Clearly  there 
is  no  need  to  settle  for  an  approximate  solution  if  an  exact  solution  is  avail¬ 
able  and  usable.  Hence  we  will  pursue  the  exact  solution  of  (86)  until  it  is 
shown  that  an  exact  solution  has  undesirable  properties  making  it  useless. 
Then  wc  will  return  to  the  second,  approximate  method. 

The  matrix  R  may  have  a  very  small  determinant  (/.  e.,  may  be  ill- 
conditioned)  in  spite  of  careful  choice  of  tjN,  so  that  its  inverse  R  '  may  be 
numerically  difificult  to  compute.  To  avoid  computational  difficulties  of  this 
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kind  and  to  define  the  p,  uniquely,  we  transform  the  trial  basis  R  to  an 
orthonormal  basis  O  by  means  of  a  Gram-Schmidt  orthogonali/ation  pro¬ 
cedure.  fhus  we  obtain  a  set  of  mutually  orthogonal  basis  vectors 
where _/  =  1, ,  M  using 

i-l  M 

UEj)  =  -  Yj  Z  ‘  =  2, ... ,  n 


H^j) 


1  he  vectors  can  be  orthogonali/cd  in  any  order,  but  the  procedure  is  de¬ 
scribed  for  ascending  /  without  loss  of  generality.  It  is  assumed  that  the 
4>{Ej,  cJiV)  are  normalized  to  unity  before  the  production  of  the  0',.  An  im¬ 
proved  basis  is  obtained  at  each  step  of  the  orthogonali/ation,  containing 
one  additional  orthogonali/cd  basis  vector  at  each  step.  The  transpo.se  of 
each  improved  basis  can  be  represented  as  the  transpose  of  the  previous  ba¬ 
sis  multiplied  by  a  simple  lower  triangular  matrix  which  has  main  diagonal 
terms  of  unity  and  no  nonzero  off-diagonal  terms,  except  for  one  (/th)  row 
which  orthogonalizes  the  next  (/th)  basis  vector,  fhus  the  transpose  of  the 
final,  fully  orthonormalized  basis  can  be  obtained  by  multiplying  the  trans¬ 
pose  of  the  trial  basis  by  n  lower  triangular  matrices,  or  equivalently  by  one 
lower  triangular  matrix  which  is  the  product  of  the  n  matrices.  (The  product 
of  lower  triangular  matrices  is  also  a  lower  triangular  matrix.)  In  matrix  no¬ 
tation,  wc  write  the  transformation  as 

0^=TR^  (91) 

where  the  lower  triangular  matrix  T  summarizes  the  Gram-Schmidt  proce¬ 
dure.  Solving  this  equation  for  the  inverse  of  matrix  R  and  substituting  into 
equation  (89)  gives 

(92) 

By  the  associative  law  for  matrices,  this  can  be  expressed  as 

k,„,,  =  {k''r){o-'r),  (93) 


A^,.,a/=(Tk)^(0-’r).  (94) 

Because  matrix  O  is  orthogonal  its  inverse  is  numerically  accessible,  fhe 
multiplication  of  the  energy  distribution  function  f  by  matrix  O  '  projects  the 
distribution  onto  the  orthonormal  basis,  which  was  transformed  from  the 
original  trial  basis  of  equilibrium  distribution  functions.  The  first  parenthe¬ 
tical  expression  likewise  transforms  the  collision  volumes  belonging  to  the 
trial  basis  into  collision  volumes  belonging  to  the  orthonormal  basis. 

Note,  however,  that  projecting  the  energy  distribution  onto  an  orthonormal 
basis  docs  not  itself  constitute  the  proposed  method  of  calculating  swarm 
mobility  or  collision  volume.  Rather,  the  essential  feature  of  the  proposed 

method  is  the  identification  of  the  K,  in  equation  (87)  (or  in  eq  (9))  as  be- 
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Kigurc  4.  Equilibrium  energy  distributions  as  trial  basis  vectors:  l:xamplcs  of  trial  basis  vectors  used  for 
projection  onto  an  equilibrium  basis.  Lquilibrium  energy  distributions  c,l N)  shown  arc  for 
i=i  (solid  line),  10  (dotted  line),  20  (dashed  line),  30  (dash -dotted  line),  40  (dash -dashed  line),  and 
SO  (second  dotted  line).  All  curves  arc  normalised  to  a  maximum  amplitude  of  unity. 


longing  to  equilibrium  energy  distributions  c,l N)  (or  (/>(/:,  A,)  in  cq  (6) 
el  seq.).  To  demonstrate  the  projection  method,  we  calculated  «  =  50  trial 
basis  functions  for  air  for 

i-l 

e,//V  =  O.OI  ,  /  =  1 , ...  ,  50  (95) 

with  a  choice  of  /c  =  1.2,  which  gives  a  range  of  values  for  c,//V  starting  with 
0  Td,  0.01  I'd,  etc.,  and  ending  with  379  Td.  'fhc  finite  mesh  equation  was 
solved  with  M  =  550  mesh  points  from  -^kT  to  about  29  eV,  to  obtain  equi¬ 
librium  energy  distributions  for  the  50  separate  values  of  t.jN,  using  the 
cross-sections  previously  discussed  for  air,  except  that  attachment  and 
avalanching  cross-sections  were  omitted  to  simplify  the  calculation.  Results 
were  saved  on  computer  disk  for  later  use  with  the  projection  method,  be¬ 
cause  of  extensive  computer  time  required  to  calculate  the  set  of  50  distrib¬ 
utions.  Equilibrium  distributions  for  i  =  I,  10,  20,  30,  40,  and  50  are  shown 
in  Figure  4. 
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Figure  5.  Orthonormal  basis  vectors  derived  from  equilibrium  energy  distributions:  First  four  orthonormal 
basis  vectors  derived  from  4)(IC,,  eJX),  i  =  1, 2,3,4.  Vectors  I  through  4  are  in  the  order  solid,  dotted, 
dashed,  and  dash -dotted.  All  curves  arc  normalized  to  a  maximum  amplitude  of  unity. 


After  these  distributions  were  calculated  as  the  trial  basis,  an  algorithm  was, 
applied  to  orthonormalize  the  set  (using  the  Gram-Schmidt  procedure  men¬ 
tioned)  in  ascending  order  of  /.  As  a  practical  matter,  the 
orthonormali/ation  was  accomplished  with  a  weight  function  proportional 
to  the  energy  mesh  spacing  to  account  for  a  noncquispaced  energy  mesh. 
The  first  four  orthonormal  basis  vectors  so  obtained  arc  shown  in  f  igure  5. 
fypical  of  orthogonal  functions,  the  number  of  sign  changes  increases  with 
the  ordinal  of  the  basis  vector. 

Because  .V/  was  chosen  greater  than  n,  the  orthogonali/ation  must  incorpo¬ 
rate  some  quantification  of  the  notion  of  “goodness  of  fit!”  which  is  not 
needed  when  Xf  =  n.  This  is  achieved  by  selecting  the  order  in  which  the  trial 
basis  vectors  arc  orthogonalized.  A  good  though  possibly  suboptimal  or¬ 
dering  is,  at  the  k  th  step,  to  orthogonalize  the  remaining  trial  basis  vector 
who.se  orthogonal  form  has  the  largest  inner  product  with  the  residual  non¬ 
equilibrium  distribution  from  the  previous  step  The  residual  nonequi¬ 

librium  distribution  at  any  step  is  the  residual  for  the  previous  step  less  its 
inner  product  with  the  orthogonal  form  of  the  trial  basis  vector  selected  for 
orthogonali/ation  at  that  step: 
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6.  INoncquilibrium  energy  distributions  for  step-function  electric  field:  I1mc-dcpcndcnt  finite  mesh 
equation  solutions  for  a  lO-Td  step-function  electric  field,  at  0  ps  (solid  line),  lO  ps  (dotted  line),  50 
ps  (dashed  line),  200  ps  (dash-dotted  line),  I  ns  (dash -dashed  line),  4  ns  (second  dotted  line),  and  20 
ns  (second  dashed  line)  after  onset  of  the  step  function.  All  curves  arc  normalised  to  unit  amplitude. 


rP  =  F[t,Ej), 


V  V 


i=\ 


As  before,  the  prime  notation  denotes  the  orthogonalized  form  of  the  basis 
vector.  Thus,  at  the  Ath  step,  we  orthogonalize  the  previously 
unorthogonalized  basis  vector  0(Zi„c,//V)  for  which  the  magnitude  of  the  in¬ 
ner  product 

l=\ 


is  maximized.  Selecting  the  order  of  orthogonalization  in  this  way  requires 
that  a  trial  orthogonal  form  be  calculated  for  each  vector  tested  using  the 
inner  product.  Because  any  orthogonal  form  depends  on  the  prior  order  of 
orthogonalization,  the  trial  orthogonal  form  may  not  be  the  same  as  the  final 
orthogonal  form  calculated  for  a  given  basis  vector. 
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I'igurc  7.  Noncquilibrium  energy  tfiMributions  for  !(quarc-wavc  electric  field:  rime-dependent  finite  mesh 
equation  solutions  for  a  100-  I  d  square-wave  electric  field,  at  0  ps  (solid  line),  10  ps  (dotted  line),  50 
ps  (dashed  line),  200  ps  (dash-dot  line),  I  ns  (dash-dashed  line),  6  ns  (second  dotted  line),  and  1 1  ns 
(second  dashed  line)  after  onset  of  the  square  wave.  Ilie  square  wave  lasts  1  ns.  All  curves  are 
normalized  to  unit  amplitude. 


We  used  this  ordering  of  the  orthogonalization  process  and  the  consequent 
order  of  rows  of  matrix  T  in  the  ensuing  calculations,  fhe  orthogonalization 
order  (row  order  of  matrix  T)  may  change  from  one  time-step  to  another, 
thereby  introducing  discontinuities  in  the  calculated  collision  volume. 

3.2  Obtaining  Time-Dependent  Energy  Distributions 

fo  illustrate  calculation  of  time-dependent  electron  energy  distributions,  we 
select  two  prescriptions  of  the  ambient  electric  field;  (1)  a  step  increase  from 
zero  field  to  a  constant  amplitude  of  10  Td,  and  (2)  a  square  wave  of  1-ns 
duration  and  lOO-'fd  amplitude.  An  air  density  one-thousandth  that  at  sea 
level  is  selected  to  reduce  the  collision  rate  and  emphasize  the  nonequilib¬ 
rium  aspect  of  the  calculation.  Using  time  leaps  of  10  ps,  we  solved  the  finite 
mesh  equation  for  air  for  both  the  prescribed  electric  fields.  Initial 
(Maxwellian)  and  subsequent  distributions  are  shown  in  I'igure  6  on  page 
32  for  the  first  ease  of  a  step  increase  in  the  field.  The  figure  show's  rapid 
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heating  of  the  electrons  on  this  short  time  scale,  relaxing  toward  an  equilib¬ 
rium  distribution  corresponding  to  a  constant  electric  field- value  of  10  Td. 

I'or  a  short  square-wave  electric  field  (shown  in  f  igure  7  on  page  33)  initial 
(Maxwellian)  and  subsequent  distributions  show  the  expected  maximum 
heating  at  I  ns  after  the  onset  of  the  field,  followed  by  very  slow  cooling 
back  to  the  zero  field  equilibrium.  As  the  distribution  cools,  the  high-energy 
tail  erodes  rapidly  between  2  and  3  eV,  where  downscatter  cross-sections  arc 
greatest  for  nitrogen.  A  small  island  of  electrons  lingers  above  3  eV,  where 
downscatter  cross-sections  arc  smaller.  A  striking  feature  of  the  distrib¬ 
utions  when  the  electric  field  has  returned  to  zero  is  the  growth  of  sharp 
maxima  and  minima  where  specific  resonances  transfer  electrons  from  one 
fixed  energy  to  a  lower  fixed  energy.  I'hcsc  sharp  peaks  and  valleys  occur 
because  the  electric  field  no  longer  dirfuscs  (“smooths  out”)  the  electron  en¬ 
ergy.  Although  clastic  scattering  dififuscs  electron  energy,  its  cflcct  is  much 
weaker  than  the  effect  of  the  electric  field,  fhe  presence  of  a  small  electric 
field  after  the  main  pulse  would  keep  electron  energies  well  mixed  and  pre¬ 
vent  such  an  irregular  distribution  function. 

3.2  Projecting  Time- Dependent  Distributions  onto  an  Equilibrium  Basis 

fhe  spectral  coefficients  p,  for  projecting  onto  the  orthogonalized 

equilibrium  basis  arc  plotted  for  case  1  (step-function  electric  field)  in 

figure  8  on  page  35  for  the  same  times  for  which  the  energy  distribution 
was  shown  in  f  igure  6.  The  spread  in  the  spectrum  gives  a  measure  of  how 
far  the  nonequilibrium  energy  spectrum  has  departed  from  an  equilibrium 
state,  fhe  spectral  projection  is  obtained  from  the  equation 

p  =  T^O"'f.  (96) 

'fhe  spectral  projection  for  case  2  (square-wave  electric  field)  is  plotted  in 
figure  9  on  page  36.  In  both  figures  the  dominance  of  a  single  equilibrium 
distribution  is  observed  for  every  nonequilibrium  distribution  calculated. 
Significant  additional  spectral  content  occurs  for  a  few  cases,  primarily  the 
lower  energy  cases  where  heating  has  just  begun.  At  higher  energies  the 
dominance  of  a  single  equilibrium  distribution  is  more  pronounced,  fhe  in¬ 
itial  (Maxwellian)  distribution  for  zero  electric  field  has  a  single  component 
at  vector  ordinal  1,  as  expected,  cofesponding  to  the  zero-field  equilibrium 
distribution.  At  10,  50,  and  200  ps,  dominant  components  at  vector  ordinals 
2  and  3  arc  seen,  with  2  to  4  nearby  components  of  amplitude  10  to  40  per¬ 
cent  as  great  as  the  dominant  component.  At  later  times  secondary  com¬ 
ponents  appear  to  contribute  less  and  less,  suggesting  that  the 
time-dependent  distribution  is  relaxing  toward  some  equilibrium  state. 

fhe  two  distributions  occurring  5  and  10  ns  after  the  end  of  the  square-wave 
electric  field  arc  seen  to  have  cooled  significantly  from  the  distribution  at  the 
end  of  the  electric  field  pulse,  with  dominant  vector  ordinals  of  32  and  31 
compared  to  40  for  the  latter.  I'hcsc  ragged  distributions  (as  shown  in  fig¬ 
ure  7)  arc  as  easily  represented  by  a  dominant  equilibrium  distribution  and 
much  smaller  secondary  components  as  arc  the  smooth  distributions.  The 
convergence  of  the  spectral  projection  appears  rapid  regardless  of  sharp 
peaks  and  valleys  in  the  distribution,  in  these  examples. 


34 


l-'igurc  8.  Spectral  projections  for  step-function  electric  field:  Spectral  projections  of  lime-dependent  finite  mesh 
equation  solutions  for  a  lU  ld  step-function  electric  field,  at  0  ps  (solid  line),  10  ps  (dotted  line),  50 
ps  (dashed  line),  200  ps  (dash-dotted  line),  I  ns  (dash-dashed  line),  4  ns  (dotted  line),  and  20  ns 
(dashed  line)  after  onset  of  the  step  function.  Abscissa  values  are  ordinal  numbers  of  the 
orthonormalized  equilibrium  basis  vectors  upon  which  the  solutions  are  projected.  A  basis  is  con¬ 
structed  for  each  solution  according  to  the  order  of  orlhogonalizalion  described  in  the  text,  so  that 
each  basis  is  dilfcrcnl  for  the  solutions  shown.  I  lowcvcr,  temperature  is  a  monolonically  increasing 
function  of  the  dominant  basis  vector.  All  curves  arc  normalized  to  unit  amplitude. 


3.3  Defining  a  Test  of  the  Equilibrium  Projection  Method 

It  is  important  to  define  carefully  what  test  can  be  applied  to  the  method  of 
projection  onto  an  equilibrium  basis  (cq  (9)),  so  that  a  meaningful  compar¬ 
ison  can  be  made  with  the  method  of  equation  (5).  If  an  exact  energy  dis¬ 
tribution  is  calculable  for  any  equilibrium  or  nonequilibrium  case,  depending 
on  the  electric  field,  then  the  use  of  equation  (5)  to  evaluate  the  collision 
volume  is  limited  only  by  the  accuracy  of  the  momentum  exchange  cross- 
section  data  used.  Likewise,  the  use  of  equation  (9)  is  limited  by  the  accu¬ 
racy  of  the  ensemble  collision  volume  data  used.  (We  neglect  error  arising 
from  discretization  of  the  energy  mesh,  which  in  principle  can  be  forced  to 
zero  if  a  sufficiently  fine  energy  mesh  is  used.)  Unfortunately,  some  exper¬ 
imental  error  is  present  in  undetermined  amounts  in  both  momentum  ex¬ 
change  data  and  ensemble  (swarm)  collision  volume  data.  Wo  perceive  no 
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r  igurc  9.  Spectral  projections  for  square-wave  electric  field:  Spectral  projections  of  time  dependent  finite  mesh 
equation  solutions  for  a  1 00- I  d  square wave  electric  field  of  I  ns  duration,  at  0  ps  (solid  line).  10  ps 
(dotted  line),  50  ps  (dashed  line),  200  ps  (dash-dotted  lii^^  I  ns  (dash-dashed  line),  6  ns  (dotted  line), 
and  I  1  ns  (dashed  line)  after  onset  of  the  square  vvavc.^Tibscissa  values  arc  ordinal  numbers  of  the 
orthonormali/cd  equilibrium  basis  vectors  upon  which  the  solutions  arc  projected.  All  curves  are 
norm;ili/cd  to  unit  amplitude. 


basis  here  to  say  which  method  is  more  accurate,  i'he  proposed  method  is 
not  necessarily  more  convenient  because  it  uses  ensemble  collision  volume 
data,  since  computation  of  the  energy  distribution  function  requires  use  of 
the  momentum  e.xchange  cross-.section  anyway. 

On  the  other  hand,  equation  (9)  offers  the  possibility  of  variational  accuracy 
in  estimates  of  ensemble  collision  volume  when  there  are  errors  in  calcu¬ 
lation  of  the  energy  distribution  function  and  in  momentum  exchange 
cross-section  data.  What  is  needed  to  make  a  comparison,  then,  is  to  (1) 
assume  that  a  baseline  set  of  cross-section  data  and  its  derived  energy  dis¬ 
tribution  functions  are  exactly  correct,  (2)  calculate  collision  volume  using 
cither  or  both  methods  (results  must  be  the  same)  from  the  baseline  cross- 
section  data  and  derived  energy  distribution  functions,  to  serve  as  a  standard 
for  comparison,  (3)  introduce  “error"  to  the  baseline  cross-section  data  and 
derive  perturbed  energy  distribution  functions,  to  serve  as  approximate  data, 
and  (4)  calculate  collision  volume  using  both  methods  from  the  approximate 
data  (results  will  be  different)  and  compare  results  with  the  standard  ob- 
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taincd  in  step  (2).  In  this  way  we  can  estimate  the  ability  of  the  proposed 
method  to  reduce  eHects  of  cross-section  error. 

I'he  cross-section  data  used  for  calculations  up  to  this  point  will  describe  the 
baseline  cross-section  standard,  and  equilibrium  energy  distribution  func¬ 
tions  obtained  using  the  finite  mesh  equation  will  be  assumed  to  be  correct 
for  the  purposes  of  the  comparison.  linsemble  collision  volume  obtained 
using  equation  (4)  from  these  cross-sections  and  energy  distribution  func¬ 
tions  will  describe  standard  equilibrium  collision  volumes,  assumed  to  be 
exactly  correct  equilibrium  values  Ibr  each  electric  field  strength  c.jX  in  the 
basis. 

We  define  two  erroneous  cross-section  sets  incorporating  din'erent  types  of 
error  to  test  the  ability  of  the  proposed  method  to  reduce  mobility  error  due 
to  error  in  the  energy  distribution  function.  In  set  A  every  cross-section  is 
twice  as  large  as  in  the  baseline  set.  This  causes  an  energy  scale  factor  error 
of  a  factor  of  two  in  the  calculated  energy  distribution  function.  In  set  B  all 
rotational  cros.s-.scctions  arc  twice  as  large  as  in  the  baseline  set.  fhis  causes 
a  more  complex  type  of  error  in  the  calculated  energy  distribution  function. 

3.4  Orthogonal  Projection  Results  for  a  Time- Dependent  Distribution 

Using  the  test  defined  in  the  previous  section,  we  calculated  collision  volume 
for  time-dependent  energy  distributions  using  cross-section  set  A.  ('ollision 
volume  was  then  obtained  using  the  orthonorrhal  projection  method  de¬ 
scribed.  I'he  orthonormal  projection  method  proved  to  be  generally  unsat¬ 
isfactory,  leading  to  large  positive  and  negative  calculated  collision  volumes 
for  the  successive  time-dependent  energy  distributions.  Scrutiny  of  the  ele¬ 
ments  of  the  orthogonali/ing  matrix  T  showed  that  the  maximum  si/c  of  el¬ 
ements  in  a  row  increased  sharply  as  the  ordinal  of  the  corresponding  basis 
vector  increased.  Thus  the  dominant  row  maximum  was  about  10'^  times 
larger  than  the  least  dominant  row  maximum,  'fhis  means  that  the  elements 
of  the  least  dominant  row  arc  very  large  and  of  varying  sign,  leading  to  very 
large  inner  products  with  collision  volume  vector  k.  fhis  is  a  recognizable 
consequence  of  the  ill-conditioning  of  the  linear  system  in  equation  (86). 
Although  the  orthogonali/ing  matrix  T  controls  error  growth  in  the  calcu¬ 
lation  of  matrix  O,  it  exacerbates  error  growth  in  the  product  Tk  in  equation 
(94).  I  his  is  due  to  the  occurrence  of  both  positive  and  negative  row  ele¬ 
ments  in  R,  T,  and  O  . 

It  was  observed  that  the  spectral  form  (eq  (87))  for  an  orthonormal 
projection  resembles  an  asymptotic  series  at  later  times  when  the  energy 
distribution  is  close  to  equilibrium.  An  asymptotic  limit  is  approached  for 

the  value  of  after  a  few  terms  in  the  expansion  (in  order  from  most 
dominant  basis  vector  to  least  dominant).  Unfortunately  this  is  not  the  case 
at  intermediate  times  when  the  distribution  is  far  from  equilibrium. 

3.5  Nonnegative  Projection  Method 

A  stabler  system  results  from  requiring  row  elements  in  R  to  he  nonnegative, 
so  that  each  p,  is  nonnegative  also,  fhe  projection  is  then  approximate,  be¬ 
cause  the  exact  solution  requires  no  restriction  on  the  sign  of  elements  of 
R.  Such  a  system  may  be  solved  by  defining  the  solution  as  the  minimum 
of  a  penalty  function  which  measures  the  weighted  error  raised  to  some 
power.  The  penalty  function  is  minimized  using  a  suitable  nonlinear  opti¬ 
mization  procedure,  yielding  the  desired  solution  for  p,. 


37 


1 


Wc  define  an  approximation  to  a  general  energy  distribution  function  /•'(/:^) 
as 


i=i 


where  the  coefficients  n,  may  be  positive  or  negative,  but  their  squares  are 
obviously  nonnegative.  Hence  the  projection  (^7)  is  nonnegativc.  I'hc  pen¬ 
alty  function  is  defined  as 


/>{7rd  = 


Wj 


-|2 


(=1 


(98) 


f  or  any  time-dependent  energy  distribution  function  /'(/:,),  minimization  of 
E  over  the  allowable  space  of  n,  yields  a  solution  for  the  n,.  I'hc  solution 
may  be  unique,  depending  on  whether  the  minimum  is  global  or  merely  local. 
Wc  multiply  the  bracketed  term  by  a  weight  function  w,  to  emphasize  sig¬ 
nificant  segments  of  the  distribution,  for  calculation  of  ensemble  collision 
volume,  an  appropriate  choice  of  weight  function  is  the  collision  volume 
raised  to  the  same  power  as  the  bracketed  quantity.  Minimizing 
such  a  penalty  function  yields  a  w'cightcd  least-squares  approximation  to 
E{E,)  and  a  least-squares  approximation  to  the  ensemble  collision  volume. 

Wc  have  obtained  such  a  nonnegativc  projection  approximation  to  the  en¬ 
ergy  distribution  function,  using  a  conjugate  gradient  minimization  proce¬ 
dure.  I'hc  projection  coefficients  obtained  were  sharply  defined  in  each  case, 
consisting  of  a  few  adjacent  nonzero  components  and  occasionally  a  smaller 
component  some  distance  away.  When  collision  volume  was  calculated  us¬ 
ing  the  spectral  form  for  the  baseline  cross-section  set  and  erroneous  cross- 
section  sets  A  and  H,  and  compared  with  results  using  the  conventional 
method  of  equation  (5),  the  spectral  form  gave  poor  results.  We  infer  that 
the  approximation  error  was  too  large  to  give  acceptable  results. 


The  calculation  of  collision  volume  is  improved  by  using  the  approximation 
error  to  calculate  a  residual  collision  volume  using  the  conventional  method, 
which  is  then  added  to  the  spectral  calculation  of  collision  volume.  Thus, 
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describes  an  approximate  projection  with  “cleanup”  of  the  residual  using  the 
method  of  equation  (5).  Although  improvement  was  noted,  the  calculated 
ensemble  collision  volume  still  compared  poorly  to  a  conventional  calcu¬ 
lation.  After  some  study,  this  was  attributed  to  poor  linear  independence 
of  adjacent  nonzero  components  in  the  nonnegative  projection.  That  is,  in 
equation  (22)  the  term  r{E,i)  is  not  small  compared  to  J{E,i)  for  the  case 
under  consideration. 


As  a  last  resort  the  nonnegative  projection  w'as  restricted  to  a  single  nonzero 
component  and  the  residual  was  “cleaned  up.”  This  single-vector  projection 
with  “cleanup”  is  equivalent  to  (1)  finding  the  equilibrium  distribution  clos¬ 
est  to  the  nonequilibrium  distribution  E(E„t)  (closest  in  the  sense  of  mini¬ 
mizing  the  least-square  error  previous  discussed),  (2)  applying  a  conventional 
collision  volume  calculation  to  the  residual  (“cleanup”),  and  (3)  adding  to 
the  latter  result  the  ensemble  collision  volume  of  the  selected  equilibrium 
distribution.  The  results  obtained  are  shown  in  I'igurc  10  on  page  39  for  the 


.38 


PERCENT  ERROR  (*10**  } 


TDIE  CS) 


10“ 


Figurc  10.  Krror  in  calculated  ensemble  collision  volume  for  step-function  K-f!cld:  Plotted  error  for  projection 
method  and  conventional  mctliod  based  on  various  cross-section  sets,  compared  to  exact  calculation 
of  ensemble  collision  volume  for  the  ease  of  a  square-wave  electric  field.  Dotted  line  is  error  for 
conventional  method  using  Set  A.  Dashed  line  is  error  for  projection  method  using  Set  A.  Dash-dot 
line  is  error  for  conventional  method  using  Set  B.  Dash-dash  line  is  error  for  projection  method 
using  Set  B.  Solid  line  is  error  for  projection  method  using  Baseline  Set  (exactly  zero). 


-Step- function  electric  field  and  in  I-igurc  1 1  on  page  40  for  the  square-wave 
electric  field.  Thc.se  rc.sults  arc  the  most  successful  obtained  by  the  equilib¬ 
rium  projection  method,  b'or  the  lO-'fd  step-function  electric  field,  the 
projection  method  sharply  reduces  error  (due  to  cross-section  sets  A  and  B) 
for  the  first  100  ps,  while  the  electron  energy  distribution  is  clo.se  to  its  initial 
equilibrium  state.  As  the  error  in  the  conventionally  calculated  collision 
volume  fortuitously  passes  through  zero,  the  projection  method  displays 
more  error,  and  then  both  methods  give  nearly  equal  error  at  late  times.  I'or 
the  lOO-'fd  square-wave  electric  field,  the  projection  method  gives  zero  error 
at  the  initial  equilibrium  but  promptly  develops  larger  error  at  the  next 
time-step  10  ps  later,  before  converging  to  nearly  the  same  error  as  the  con¬ 
ventional  method  at  late  times. 
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I'igurc  1 1 .  Krror  in  calculated  ensemble  collision  volume  for  square-wave  K-field:  Plollcd  error  for  projection 
method  and  conventional  method  based  on  various  cross-section  sets,  compared  to  exact  calculation 
of  ensemble  collision  volume  for  the  ease  of  a  step-function  electric  field.  Dotted  line  is  error  for 
conventional  method  using  Set  A.  Dashed  line  is  error  for  projection  method  using  Set  A.  Dash-dot 
line  is  error  for  conventional  method  using  Set  B.  Dash-dash  line  is  error  for  projection  method 
using  Set  B.  Solid  line  is  error  for  projection  method  using  Baseline  Set  (exactly  zero). 


4  Discussion 


The  use  of  an  orthogonal  basis  for  the  projection  method  was  found  to  cre¬ 
ate  unusable  results  for  the  spectral  form  of  K,„u,  because  of  the  occurrence 
of  very  large  positive  and  negative  elements  of  matrix  T.  This  is  because  of 
the  inherent  linear  dependence  of  the  basis  (^(/:„  c,//V),  and  is  a  fundamental 
shortcoming  of  the  projection  method.  If  the  basis  could  be  restricted  to  a 
small  number  of  basis  vectors  possessing  a  high  degree  of  linear  independ¬ 
ence,  the  projection  method  might  make  a  better  showing  in  cases  not  dis¬ 
cussed.  l-'igure  4  shows  that  the  equilibrium  distributions  for  /  =  1  and 
i  =  50  arc  the  most  nearly  linearly  independent  of  those  shown,  as  would  be 
expected.  However,  it  is  dilTicult  to  find  a  realistic  case  which  consists  of  a 
sum  of  these  two  distributions  alone.  Comparison  of  the  equilibrium  dis¬ 
tributions  of  f  igure  4  and  the  nonequilibrium  distributions  of  f  igure  6  shows 
that  their  shapes  compare  well  up  to  the  maximum  but  fall  off  differently  at 
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higher  energies.  The  equilibrium  distributions  have  a  smaller  high-energy 
“tail”  due  to  more  extensive  downscatter  from  rotational  and  vibrational 
excitations.  This  explains,  to  some  extent,  why  single-vector  projection  with 
“cleanup”  performs  well  on  nonequilibrium  distributions  whose  maxima  are 
concentrated  at  or  below  about  0. 1  eV.  The  poor  similarity  observed  be¬ 
tween  equilibrium  and  nonequilibrium  distributions  approaching  1  eV  at 
maximum  tends  to  prohibit  stationary  behavior  of  projection  coefficients  at 
these  energies.  In  fact,  collision  volume  error  at  these  energies  is  nearly  the 
same  for  the  projection  method  (with  “cleanup”)  and  the  conventional 
method. 

fhe  method  of  selection  of  a  dominant  basis  vector  was  the  same  for  an 
orthogonal  basis  and  for  a  nonnegative  projection  basis,  fhe  same  compo¬ 
nent  is  identified  in  cither  case.  No  assumption  was  made  about  the  form 
of  error  induced  in  the  energy  distribution  function  by  the  erroneous  cross- 
.section  set  used.  Several  other  approaches  were  tested,  including  best  fit  to 
first  and  second  energy  moments  (average  energy  and  energy  spread),  as 
ways  of  selecting  the  dominant  basis  vector.  No  significant  improvement 
was  found,  except  that  it  became  clear  that  knowledge  of  the  type  of  error 
could  be  exploited  to  produce  a  superior  error  reduction.  Deliberately  using 
a  “hotter”  basis  vector  than  indicated  caused  a  better  reduction  in  collision 
volume  error,  because  of  cancellations  in  error  contributed  by  different  parts 
of  the  energy  distribution.  This  method  was  not  our  objective,  yet  may  be 
useful  in  certain  circumstances.  Thus  if  models  of  cross-section  error  are 
defined,  the  method  of  dominant  basis  vector  selection  can  be  biased  to  re¬ 
duce  error  further  than  a  “blind”  method,  such  as  we  have  used. 

Both  orthogonal  and  nonnegative  projection  methods  work  about  equally 
well  for  energy  distributions  near  equilibrium,  so  that  significant  error  re¬ 
duction  is  possible  when  “cleanup”  is  added  (eq  (99))  to  the  spectral  form  for 
ensemble  collision  volume.  Unfortunately  both  methods  serve  poorly  when 
the  energy  distribution  is  far  from  equilibrium.  For  the  strong  100-Td  elec¬ 
tric  field,  this  state  is  attained  within  10  ps,  and  within  100  ps  for  a  10-Td 
electric  field.  I  his  may  be  due  to  the  large  rotational  and  vibrational  re¬ 
sponses  above  0. 1  cV  for  air  which,  in  equilibrium,  depopulate  that  portion 
of  the  energy  spectrum  filled  by  strong  electric  fields.  This  hypothesis  is 
strengthened  by  the  fact  that  ensemble  collision  volume  estimates  by  the 
projection  method  are  too  low  in  this  regime,  presumably  because  higher 
energies  are  underpopulated  in  equilibrium  compared  to  nonequilibrium. 
(Higher  energies  carry  a  larger  collision  volume  per  electron.)  Gases  that 
lack  such  a  strong  depopulating  process  may  yield  better  projection  method 
results,  although  such  a  case  lacks  interest.  Likewise,  distributions  substan¬ 
tially  above  2  cV  in  air  might  avoid  this  underpopulation  effect  and  yield 
better  nonequilibrium  results  for  the  projection  method. 

The  single-vector  projection  with  “cleanup”  could  have  been  derived  inde¬ 
pendently  of  the  projection-related  theory  of  this  effort,  as  it  is  conceptually 
very  simple.  It  simply  trades  most  of  the  collision  volume  calculation  using 
equation  (5)  for  an  experimentally  measured  collision  volume  represented 
by  equation  (9).  In  the  context  of  its  derivation  here,  the  variational  prop¬ 
erties  of  the  single- vector  projection  with  “cleanup”  become  apparent.  If  the 
error  associated  with  the  measurement  is  as  great  or  greater  than  the  error 
associated  with  the  energy  distribution  function,  then  nothing  is  gained  in 
accuracy.  Generally  we  expect  more  accuracy  from  an  equilibrium  meas¬ 
urement  than  from  a  nonequilibrium  calculation. 
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Although  this  dfort  is  directed  at  estimates  of  ensemble  collision  volume,  the 
same  computational  machinery  can  be  applied  to  other  electron  swarm 
properties  as  well,  such  as  the  cocfTicient  of  attachment  or  avalanching.  In 
such  a  case  a  weight  function  equal  to  that  energy-dependent  coefficient 
would  be  substituted  for  K{Ej).  fhe  previous  conjecture  regarding  under¬ 
population  at  several  electron  volts  applies  as  well  in  this  case,  where  poor 
projection  method  results  would  be  expected  until  large  spectrum  content 
develops  at  higher  energies.  Avalanching,  however,  does  not  become  im¬ 
portant  until  the  spectrum  content  develops  at  higher  energies. 

Nothing  has  been  specified  about  the  source  of  electrons  described  by  the 
energy  distribution  function,  up  to  this  point.  In  nonattaching  gas,  electrons 
will  linger  until  absorbed  by  waifs,  for  example,  f  or  attaching  ga.scs  such 
as  air,  a  continual  source  of  electrons  is  required  to  sustain  a  significant 
electron  density  for  times  much  longer  than  the  attachment  time.  While  an 
experimental  apparatus  might  generate  clouds  of  near-thermal  electrons, 
there  is  also  the  case  of  electrons  produced  by  ionizing  radiation.  Such 
electrons  arc  produced  throughout  an  exposed  volume  at  relatively  high  en¬ 
ergies  from  10  to  100  eV.  The  projection  method  highlights  the  possible 
value  of  equilibrium  measurements  of  swarm  parameters  in  a  volume  ex¬ 
posed  to  constant  ionizing  radiation.  In  such  a  measurement  the  electric 
field  would  be  imposed  on  an  ensemble  of  electrons  whose  energy  distrib¬ 
ution  incorporates  the  peculiar  spread  of  electrons  throughout  higher  ener¬ 
gies  derived  from  their  birth  through  ionizing  radiation.  Such  equilibrium 
measurements  could  shed  considerable  light  on  swaiin  properties  of  near- 
equilibrium  and  nonequilibrium  distributions  in  ionizing  radiation. 


5  Conclusions 


A  dilferential  equation  has  been  derived  describing  the  time  evolution  of  the 
energy  distribution  function  for  free  electrons  in  a  gas  in  a  transient  electric 
field.  A  finiic-dilfcrcncc  approximation  and  a  time-stepping  algorithm  using 
matrix-vector  techniques  were  devised  to  solve  the  differential  equation.  An 
extensive  published  cross-section  set  was  added  to  the  time-stepping  algo¬ 
rithm  to  enable  calculation  of  realistic  energy  distribution  functions  in  air. 
A  projection  method  was  described  which  allows  computation  of  electron 
swarm  properties  such  as  collision  volume  (related  to  mobility)  using  meas¬ 
ured  data  for  equilibrium  distribution.s,  and  the  projection  method  was 
shown  to  have  desirable  variational  properties  in  calculating  the  swarm 
properties. 

Several  projective  methods  were  compared  with  conventional  techniques  for 
their  ability  to  reduce  error  in  estimates  of  collision  volume  arising  from 
cross-section-induced  error  in  the  energy  distribution  function,  including 
orthogonal  and  nonnegative  projections.  A  single-vector  projection  with 
“cleanup”  of  approximation  error  was  found  to  perform  best,  leading  to  a 
significant  reduction  of  collision  volume  error  when  the  distribution  function 
maximum  was  not  greater  than  about  0. 1  eV.  Other  projections  performed 
poorly  compared  to  conventional  methods  of  calculating  swarm  collision 
volume,  because  of  the  inherent  linear  dependence  of  the  projection  basis  of 
equilibrium  distributions  and  the  resulting  ill-conditioning  of  the  projection 
matrix.  If  a  model  of  the  energy  dependence  of  the  cross-section  error  is 
known,  it  may  be  possible  to  exploit  that  knowledge  to  bias  the  single-vector 
projection  to  give  better  results. 
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6  Recommendations 


I'lcctron  mobility  estimates  for  near-thermal  (with  a  peak  not  greater  than 
0.1  cV)  nonequilibrium  energy  distributions  can  be  made  with  greater  accu¬ 
racy  using  a  single-vector  projection  (onto  an  equilibrium  basis)  with 
“cleanup.”  Estimates  based  on  the  mobility  of  an  equilibrium  distribution 
having  the  same  average  energy  as  the  nonequilibrium  distribution  arc  a 
special  case  of  the  above  without  “cleanup.” 

for  electron  .swarms  populated  by  ionizing  radiation,  nicasurcmcnfs  of 
equilibrium  swarm  properties  in  a  time-constant  electric  field  and  radiation 
source  would  likely  be  very  useful  and  shed  new  light  on  swarm  properties 
of  nonequilibrium  electron  distributions  arising  in  this  way. 
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